{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Different objective functions for Approximate GPs\n",
    "\n",
    "\n",
    "(*N.B.* this tutorial assumes that you are familiar with **inducing point methods**.\n",
    "For an introduction to these methods, please see [Quinonero-Candela and Rasmussen, 2005](http://www.jmlr.org/papers/volume6/quinonero-candela05a/quinonero-candela05a.pdf).)\n",
    "\n",
    "## Overview\n",
    "\n",
    "Approximate Gaussian processes learn an approximate posterior distribution\n",
    "\n",
    "$$ p(\\mathbf f(X) \\mid y) \\approx q(\\mathbf f) = \\int( p(\\mathbf f \\mid \\mathbf u) \\: q(\\mathbf u) d \\mathbf u$$\n",
    "\n",
    "where $\\mathbf u$ are the function values at some **inducing points** $\\mathbf u$.\n",
    "Typically, $q(\\mathbf u)$ is chosen to be a multivariate Gaussian -- i.e. $q(\\mathbf u) = \\mathcal N (\\mathbf m, \\mathbf S)$.\n",
    "\n",
    "We choose the approximate posterior $\\int( p(\\mathbf f \\mid \\mathbf u) \\: q(\\mathbf u) d \\mathbf u$ by optimizing the parameters of $q(\\mathbf u)$ (i.e. $\\mathbf m$ and $\\mathbf S$. There are several objectives that we can use for optimization. We'll test out the following two:\n",
    "\n",
    "1. The Variational ELBO (see [Hensman et al., 2015](http://proceedings.mlr.press/v38/hensman15.pdf))\n",
    "2. The Predictive Cross Entropy (see [Jankowiak et al., 2019](http://arxiv.org/))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experimental setup\n",
    "\n",
    "We're going to train an approximate GP on a 1D regression dataset with heteroskedastic noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Text(0, 0.5, 'y'), Text(0.5, 0, 'x')]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAADNCAYAAAArFbJhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAVpElEQVR4nO3dfWzcd2HH8ffP13MYD2VphAdoLMPsQVtmiuQv2nmUdA3s4dJpUulUKXG38LBuQ4MFqsp0mtuQYiB4gRSQNlBTL5DDGRVsKhu2EcJVvVg+xNdMmxlimzapi9hDplCrMIjvbP/2h/07zpev7353vvs93H1eUqXYPt/v++3dffx9/nm+7yMiIjv1xV0AEZEkUjiKiDgoHEVEHBSOIiIOCkcREQeFo4iIw01xFyCMBx98UOuNRKQjzpw547m+n4pwBDh9+nTTv3P16lUGBgY6UJpoqR7J0y116fV6nDp1atefqVstIuKgcBQRcVA4iog4pGbMUUSkWrFYZGFhgcOHDzM4ONj251c4ikjqFItF8vk8pVKJ/v5+pqenyefzbb2GutUikjoLCwuUSiU2NjYolUosLS21/RoKRxFJncOHD9Pf308mk6G/v5+RkZG2X0PdahFJnVwux+zsrMYcRURq5XI5crkcsLUIvN3UrRYRcVA4iog4KBxFRBwUjiIiDgpHEREHhaOIiIPCUUTEQeEoIuKgcBQRcVA4iog4KBxFRBwUjiIiDgpHEREHhaOIiIPCUUTEQeEoIuKgcBSR1CgWi0xOTlIsFjt+rVhOAjfGvAqYAL4O/DhwzVr7SBxlEZHmVd8WNTiNO4prVt9xcHZ2tqPXjus2CbcAf2mtfRLAGPNNY8wXrbXLMZVHREKKOqQCtXccXFhY6L5wtNZ+reZbfcD/xVEWEWlO1CEVCO44GITy4cOHO3q92G+wZYy5C/iStfZb9R7Xyg10VldXWy1WoqgeydMtdWmlHkNDQ2SzWQCy2SxDQ0MducFVrcHBQaanp1laWmJkZITBwcHKdTvxesQajsaYO4A7gHc1euzAwEBL12j195JG9UiebqlLs/XI5/PMzc1FPuYYXDufzzt/1u7XI7ZwNMbcCbweOAm8zBhz0Fq7FFd5RCS86tuidqtYlvIYY4aBzwI54CngSeBn4yiLiIhLXBMyy8AL47i2iEgYWgQuIuKgcBQRcVA4iog4KBxFRBwUjiLSMVEeFNFuse+QEZHu1Moe7DgOtNiNwlFEOqLZPdhxHWixG3WrRaQjgoMiMplMqIMiXGEaJ7UcRaQjcrkcs7OzobvJUZ+604jCUUQ6ppk92M2GaacpHEUkMZJ0oIXGHEVEHBSOIiIOCkcREQeFo4iIg8JRRMRB4Sgi4qBwFBFxUDiKiDgoHEVEHBSOIhK7JJ77qO2DIhKrpB1VFlDLUURilbSjygIKRxGJVbPnPkZF3WoRiVXSjioLKBxFJHZJOqosoG61iIiDwlFExEHhKCLiENuYozHmpcAEcKu19rVxlUNExCXOCZnbgCeB13TiyYvFIjMzMxw9ehQgcTNhIpJssYWjtfZzxphf7sRzV6+4f/TRR/E8j/X1dTKZDCdOnGB0dBRQYIo0UiwWe/ZzkpqlPFevXg392JmZmcqK+83NTQB832djY4Pz58/zqU99CoCNjQ0ymQz33HMPd999NwBLS0uMjIwwPDzc/kq0aHV1Ne4itEW31AO6py716rG8vMyxY8col8tks1kuXbqUqM9FtU68HqkJx4GBgdCPPXr0KB//+McplUpkMhk8z6NUKuH7Pr7vUy6XgR8G5mc+8xmeeOIJZwszKX8tm6l/knVLPaB76rJbPVZWViiXy2xsbFS+zufzURatKe1+PVITjs0IVtxXjzkWCgUuXrxYCb9GgXn+/HkKhQJnz57l2rVrPdmtkN4WbOsLDoRIyra+qMQ5W3078NvAy4wx48CHrbU/aNfz53I5BgcHK39Ncrkc9957b2X8BBoH5traGidPnsT3/US2JkU6Kanb+qIS54TM08DTUV6zdotSo8D0PI/NzU02NzfVmpSelMRtfVHpym51M+oF5oEDB3jggQe4fv26szWZpLPnRDqt12auez4cXaoD89ChQ7u2JtfW1piYmGB8fLwn3izSu5J6IG0nKRwbCIKytjW5trbG5uYm8/PzLC4uqqstXc11IG2r7/O0tEAVjiHVtiYnJiaYn5+vtCDV1ZY0CQJqaGgo1PKcds1cp6kFqnBsQS6XY3x8nMXFRUql0o6udqlUolAopOIvo/Sm6oDKZrPMzc01fJ+2a+a6nS3QTlM4tqj6zRJ0tYNF59XLg7T8R5KmOqCCr8O8PxvNXIfpLqdp7WTDcPQ87zd93/9CFIVJm9qu9sLCAleuXGFqaoqNjY0dy3+S3H2Q3lIdUNls1hlQzY4LNuouVz9fWtZOhmk5jnuedxvwad/3v9HpAqVVEJTFYpFCobBj+Y+62pIk1b2eoaGhG96PrYwL1usuu55vbGysY/VrlzDheAJ4Bniz53lvB+Z83/+bzhYrvYI3Xu3um+DfSR+Elt4Q/DF3HejSyrhgve5ymsYZq4UJxwywAawBvwQc9DzvV4G/833/iU4WLq1ql/9Ud7W1NlKSrpVxwXoTNmkaZ6wWJhwLwIuAvwXu8X3/XwE8z3sEUDjWUdvVrl0bqRakJFGrM9O7TdikdY92mHD8Z+B3fd//bvANz/P6gf0dK1WXCd4c1WsjNQ4pSdbuPdVp3KMdJhyP+76/Uf0N3/dLwDs7U6TuVLs2UuOQ0k5p2XWSJg3DsTYYpXXV3Yvqccg0DVJL8qRp10ma6NasEcvlcoyNjTE6Okp/fz+ZTIZMJsOVK1coFotxF09SyDUbLHuncIxJ0Ip8y1vegud5TE1Nkc/nefzxx5mcnFRQSmjBbHAmk0nVbHDSaftgjHK5HAsLC6yvr1eW+egAC2lWWmeDk07hGLPqNWC1B1hoHFLCSuNscNIpHGNW7wCLYBxycHAw7mKK9ByFYwLsdvL41NQUhUKBhx9+mHK5rC6TSIQUjgnjGod86KGHNA4pPSuuNZwKxwTSOKTIljjXcGopTwIF45CnTp3iox/9qJZpSM+Kcw2nWo4JVT0O+fKXv5yVlRUOHDhQeXOo9Si9IM4TfRSOKTA8PMz+/ft3dC90t0PpBXGu4VQ4pkR190KLxaWXxLWGU2OOKVG9Rayvr4/Nzc0dh+dqu2F6FIvFSLeIRn29bqGWY0q4Fovr8Nz0iXr2VSf2tE4txxQJTvR529vexuzsLEeOHKm0InUaSzpEPfuqE3tap3BMqeDw3H379lWW+Rw4cEDdp4SL+gQdndjTuti61caYNwJvAq4CvrX2dFxlSavd9mWr+5RcUc++6sSe1sXScjTGPB/4BPBua+17gVcbY94QR1nSLuhqX7t2bcdstiZpkit4zaIKqrDX08TNTnG1HEeAZ6y1a9tfLwJ3Al+JqTypF3SfNEkjrdDEzY3iCscB4LtVXz+3/b1duW4+3sjq6mrTv5NEYeoxODjI9PQ0586d4/Lly5VJmpmZmcQcedYtrwekoy7Ly8ssLS0xMjLC8PCw8zGrq6ssLy9z7ty5yh/WpL1vGlleXuapp57ijjvu2LWerYgrHK+ydS/swM3b39vVwEDd7Gz77yVNmHrk8/kbdtIcPHiQCxcuJGa8qVteD+hcXdpxCk2xWOT48eMNW4LLy8scP368Eox9fX309/dz9OjRVLxW1fV87LHH2trijSscl4CDxph9213r1wF/FlNZuoomadKtXd1b1xIe1/MsLS1RKpUqwXjkyBHGx8dT8z4JW89WxDIhY639PvB24GPGmAngH621Gm9sE9ckTalUolAoaMA94dq1LjHsEp6RkZHK4/bt29d0MMY9idPJpUqxLeWx1n4Z+HJc1+8F1SeaZDIZLl68yPr6ulqRCRbmFJp63e7qn4VZwjM8PNzyUp9OT+KEGV4IekozMzMcPXq0rdfX9sEuVt3FvnLlClNTUzuW+qSp+9QrGq1LrBdIrp+NjY2FumYnu+6taCZ4c7kcg4ODbR8j1Q6ZLhd0sUdHR+nv769sN5yfnyefz6uLHTNXt7TeusR63e4wXfJ2doM72aVNwrZHtRx7RNAimZiYYH5+vrJko1AoaPdETFrpltbrdjfqkruut5flOru1csPOttd7XJyH3AYUjj0k2I+9uLioccgEaKVbWq/b3ahL7rreXtcy1nbJwwZ+o8clYdujwrHH7DYOqZt3Ra/V1lG9McJ6P4uiNRY28MM8Lq5DbgMKxx4UvOmKxSKFQqHSirxy5QrFYlEBGZEkHELRys6zesIGcBK6zY0oHHtY8GEpFApcvHiRqakpCoWCutcRirp11OnrhQ38JHSbG1E49rhcLsfCwgLr6+vqXktbhA3guLvNjWgpj+hAVBEHtRwlFV2cbtCOAyUkOgpHAXZ2cYIP8YEDB3Rv7DbReYnpo3CUHYIPcfURVvv27ePs2bMKyj3o5FY76QyFo+wQfIg3NzcB2NzcZG1tjZMnT+L7vlo9LWp16Uq9Vry66Z2lcJQdam+30NfXV9mPXX0LWH0YmxPmQAnXNjxXK352dhZA3fQOUzjKDrWH5V67du2GQ3M1m92a3Zau7DYe6WrFVx/CoG56Zykc5QauD/GhQ4cqgRl8OPVhbI/dxiNdrfjqP05J32GSdgpHCSUIwuoWTq9M0nR6bG+38UhXK766DFp+1VkKRwmtuoXTK5M0YZbghD2Zu5WtdK0eMiF7p3CU0KpbOJ7n9cQkTW2XNzj/cmhoqHJYcDMncyvo0kPhKKHVdvOCSZpuOtGntqW32314stksc3Nzddcvam1juikcpSnVLZxDhw41daLP8vIyKysriR0j262l5zr/EqiEaKsnc0uyKRylZc2c6FMsFjl27BjlcjmxY5S7tfRc519ms9lKyFe3pqtn8rVnPd0UjrInu7WOqrunABMTE5U1e0ntYjZq6VWH3dDQUKX8rpn8IPw1lpheCkfZk9rWEcA73vGOythcJpPB8zzK5bJzrV6ShGnpBWFXe4K2xhe7j8JR9qy625nP57l+/Tq+7wNUdnf4vk9fXx9HjhxJ9P2yW23paXyx+ygcpW2C1lMQjJ7nkc1m8TyvMsNbLxjbvdi60fO183quFvTk5KTGGlNM4ShtU7vs5cSJE4yOjgLcME5Xq93nHTZ6vk6cr1jbgtahEOmmcJS2abTTo96d7lods3Md6QVbE0DBnmTX83VyjFDjj91B4ShtFeWYnetIr5tuuinUBFAnxwg1/tgdFI4Si9rxvjAzxbW/4zrSq1wuA40ngFq5Xlha39gdIg9HY0wfcB/wPuCItfYbUZdB4lU7Jld9us/Y2FjlMdXh4hrHcx3pFbQc19fX6e/vrzsBVK+Vu9dxQ61vTL84Wo63Al8Fvh/DtSUBGp3uAzcuqHaN442NjTmP9AqusZdWm8YNJfJwtNb+PYAxJupLS0I0Ot0H2BGeExMT3HXXXbueebjbYu1qzd5RUeOG0pFwNMZ8Cfgxx48ettZ+oZXnrDfTuZvV1dVWLpU43VaPwcFBpqenWVpaYv/+/Zw+fZpyuUw2m2VoaAiAbDaL7/tsbm4yPz/P5cuXOXXqFM8++ywjIyMMDg6Gfk8sLy9z7NixyvhkMElz6dIlhoeHnb9TXUbX9Z5++mlWVlYYGRnZ9TnSoNveW+3UkXC01v5au59zYGAg0t9Lmm6rRz6fJ5/PAzAyMnJDN3hubo6JiQnm5+crEy3lcplHHnmEYrHIhQsXQnebV1ZWKrPX8MOJm5WVlUoZXKrLWK1YLHLfffcl+hCNZnTbe6tdNFstsXN1jXO5HOPj4ywuLu7o2rYyUdLoXizNWlhYoFwuazyyy8UxW70f+EPgxcDvGWOmrbXFqMshyedaEjM5Odn0REmje7E06/Dhw2SzWQCNR3axOCZkngUmtv8Tqau2VdnqREk7l9bkcjkuXbqU6IN7Ze/UrZZUScoC6+Hh4brjlZJ+CkdJHS2wlij0xV0AEZEkUjiKiDgoHEXYWrs4OTlJsaiFE7JFY47S83Q4rbio5Sg9z3XIhIjCUXpesHYyk8loUbdUqFstPS8payclWRSOImjtpNxI3WoREQeFo4iIg8JRRMRB4Sgi4pCaCZlTp07FXQQR6SGe7/txl0FEJHHUrRYRcVA4iog4KBxFRBwUjiIiDqmZra7HGPNG4E3AVcC31p6u+fnzgLPAt4GfBs5Ya/8l8oI2EKIe7wFeCvw3MAw8bK39VuQFbaBRPaoeNwoUgBdZa78XYRFDCfF6eMA7t7/8SeBHrbVvjbSQIYWoyyvZ+ox8DXgNMG2t/ULkBa3DGPNStm7Md6u19rWOn/cBHwC+BxwEHt/LnU1T33I0xjwf+ATwbmvte4FXG2PeUPOwdwH/Ya39IHAOeDzaUjYWsh4vBO631n4I+Dzwp9GWsrGQ9cAY83PAz0dcvNBC1uNeYNVa+zFr7f3AoxEXM5SQdRkDLltrzwAfAj4cbSlDuQ14EvB2+fk9wM3W2gngPcCnjTGZVi+W+nAERoBnrLVr218vAnfWPOZOYAnAWrsC3GqMuTm6IobSsB7W2oestcHaqz62/kImTcN6bH9YxwBnizIhwryvRoFbjDF/ZIwJWixJFKYu/wO8ZPvfLwGWIypbaNbazwHfrfOQ6s/5d4DrwKFWr9cN4TjAzv9hz21/r9nHxC10GY0x/cAJYDyCcjUrTD3eD7zPWluKrFTNC1OPg2y1VD4GXADm9tJS6aAwdfkI8IvGmI8ADwN/EVHZ2qmtn/NuCMerwIuqvr55+3vNPiZuocq4HYx/DvyJtfbfIipbM+rWwxjzCmA/cI8x5sHtb99vjDHRFTGUMK/Hc8BXAbbHsG8GXhFJ6ZoTpi4XgPPbwwN3AZ81xtwSTfHapq2f824IxyXgoDFm3/bXrwO+aIy5parr/EW2uhYYY4aAf7DWPhd9UetqWA9jzI8AnwQ+Yq1dNsbcHVNZ66lbD2vtFWvtm621Z7bHt2CrPjae4u4qzPvqK8AgwPb3MmxNliVNmLq8Aviv7X8/C2ySgnwwxrzAGBMMB1R/zm8Bngf8U6vP3RXbB40xvwL8FvC/QNlae9oYMwl8x1p7ZjtUzrL14v8U8IGEzlY3qsdfAb8A/Of2r7zANWsXt0b12H7MS4DfB963/d8nrbXfjqvMLiFejxcDk8AzwKuAz1trZ+Ir8e5C1OU2tiYuvw68Eli21n4ivhLfyBhzO/A7wK+z1Xv6MPBWYMha+wfbs9UfBL4P/ATw2F5mq7siHEVE2i3xzWYRkTgoHEVEHBSOIiIOCkcREQeFo4iIg8JRRMRB4Sgi4qBwlK7ged4fe573A8/zbvc8737P82Y8z/uZuMsl6aVF4NI1PM97iK2dESXgAd/3fxBzkSTF1HKUbvJ+4PXANxWMsldqOUrX8DzvbrZOZRkDfsP3/X+PuUiSYmo5SlfwPO+twIPALFsHtf6153m3x1sqSTO1HEVEHNRyFBFxUDiKiDgoHEVEHBSOIiIOCkcREQeFo4iIg8JRRMRB4Sgi4vD/bbswq0TxgVMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define some training data\n",
    "train_x = torch.linspace(0, 1, 100)\n",
    "train_y = torch.cos(train_x * 2 * math.pi) + torch.randn(100).mul(train_x.pow(3) * 1.)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(5, 3))\n",
    "ax.scatter(train_x, train_y, c='k', marker='.', label=\"Data\")\n",
    "ax.set(xlabel=\"x\", ylabel=\"y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a simple approximate GP model, and a script to optimize/test the model with different objective functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ApproximateGPModel(gpytorch.models.ApproximateGP):\n",
    "    def __init__(self, inducing_points):\n",
    "        variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(inducing_points.size(-1))\n",
    "        variational_strategy = gpytorch.variational.VariationalStrategy(\n",
    "            self, inducing_points, variational_distribution, learn_inducing_locations=True\n",
    "        )\n",
    "        super().__init__(variational_strategy)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "training_iterations = 2 if smoke_test else 50\n",
    "\n",
    "\n",
    "# Our testing script takes in a GPyTorch MLL (objective function) class\n",
    "# and then trains/tests an approximate GP with it on the supplied dataset\n",
    "\n",
    "def train_and_test_approximate_gp(objective_function_cls):\n",
    "    model = ApproximateGPModel(torch.linspace(0, 1))\n",
    "    likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "    objective_function = objective_function_cls(likelihood, model, num_data=train_y.numel())\n",
    "    optimizer = torch.optim.Adam(list(model.parameters()) + list(likelihood.parameters()), lr=0.1)\n",
    "\n",
    "    # Train\n",
    "    model.train()\n",
    "    likelihood.train()\n",
    "    for _ in range(training_iterations):\n",
    "        output = model(train_x)\n",
    "        loss = -objective_function(output, train_y)\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        optimizer.zero_grad()\n",
    "\n",
    "    # Test\n",
    "    model.eval()\n",
    "    likelihood.eval()\n",
    "    with torch.no_grad():\n",
    "        f_dist = model(train_x)\n",
    "        mean = f_dist.mean\n",
    "        f_lower, f_upper = f_dist.confidence_region()\n",
    "        y_dist = likelihood(f_dist)\n",
    "        y_lower, y_upper = y_dist.confidence_region()\n",
    "\n",
    "    # Plot model\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(5, 3))\n",
    "    line, = ax.plot(train_x, mean, \"blue\")\n",
    "    ax.fill_between(train_x, f_lower, f_upper, color=line.get_color(), alpha=0.3, label=\"q(f)\")\n",
    "    ax.fill_between(train_x, y_lower, y_upper, color=line.get_color(), alpha=0.1, label=\"p(y)\")\n",
    "    ax.scatter(train_x, train_y, c='k', marker='.', label=\"Data\")\n",
    "    ax.legend(loc=\"best\")\n",
    "    ax.set(xlabel=\"x\", ylabel=\"y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objective Funtion 1) The Variational ELBO\n",
    "\n",
    "The **variational evidence lower bound** - or ELBO - is the most common objective function. It can be derived as an lower bound on the likelihood $p(\\mathbf y \\! \\mid \\! \\mathbf X)$:\n",
    "\n",
    "$$  \\mathcal{L}_\\text{ELBO} =  \\sum_{i=1}^N \\mathbb{E}_{q( \\mathbf u)} \\left[ \\mathbb{E}_{f( \\mathbf f \\mid \\mathbf u)} \\left[    \\log p( y_i \\! \\mid \\! f_i)  \\right] \\right] - \\beta \\: \\text{KL} \\left[ q( \\mathbf u) \\Vert p( \\mathbf u) \\right]$$\n",
    "\n",
    "where $N$ is the number of datapoints and $p(\\mathbf u)$ is the prior distribution for the inducing function\n",
    "values. For more information on this derivation, see [Hensman et al., 2015](http://proceedings.mlr.press/v38/hensman15.pdf).\n",
    "\n",
    "### How to use the variational ELBO in GPyTorch\n",
    "\n",
    "In GPyTorch, this objective function is available as [gpytorch.mlls.VariationalELBO](http://gpytorch.readthedocs.io/en/latest/marginal_log_likelihoods.html#variationalelbo)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAADNCAYAAAArFbJhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deXxU5b3/3+fMkp0shCTsu6BsCqMMW1Bww6VetWprrVbt7WZbb2tv7W211qXVurXalmuv1Vqlan+1WLFsVpGELcghKLgAGlkDIYGQhCQkM3PO8/vjyckkYbLPlvC8X695QWYmM8/JnPOZ7/NdNSEECoVCoWiNHusFKBQKRTyixFGhUChCoMRRoVAoQqDEUaFQKEKgxFGhUChCoMRRoVAoQuCM9QK6wk9+8hOVb6RQKCLCI488ooW6v0+II8D999/f7d8pLy8nJycnAquJLv3lOKD/HIs6jviip8dx3333tfuY2lYrFApFCJQ4KhQKRQj6zLa6JYFAgLKyMhobG+mo/NGyLGpqaqK4sq6haRoJCQnk5eXhdPbJj0Ch6Pf0ySuzrKyMlJQUhg4diqaF9KUC4Pf7cblcUVxZ1xBCUFVVRVlZGcOGDYv1chQKRQj65La6sbGRjIyMDoUxntE0jYyMDBobG2O9FIWiT1JUVMSjjz5KUVFRxN6jT1qOQoiYC6O9nW+5Dsuy0PWufd9omtahS0ChUISmqKiIRYsW4fP5cLvdrFy5kjFjxoT9ffqk5Rhr/H4/9957L6ZpcuDAAW6++WbuuOMOHn74YYqLi2O9PIWiX1NYWIjP58M0TXw+H4WFhRF5nz5pObbl9dc1KipOtSRN04nD0bn+DxokuPrqrltxTz/9NJdeeilOp5N169axYMECbr75ZgCuvfZaXn/99a4vXqFQdIv8/Hzcbnez5Zifnx+R9+kX4lhRoTF06Kn3BwKCrgSDS0s1ILQ4VldX853vfIfx48eze/du/H4/+/fv584776S2trZZCPPy8rj00ksZOHAgH374IZMnT+7FESkUivbwer2sXLmSwsJC8vPz8Xq9lJeXh/19+oU4RpLnn3+ehQsXctttt/Hggw8ycuRIHnvsMZxOJ6mpqVx55ZUAXHrppQDk5uZSUlKixFGhiCBerxev1xvR91A+x04oKSlhxIgRAAwfPhwA0zTbfb7T6cTn80VlbQqFInIoceyEsWPHsm/fPgAOHjwIQEpKSrvPr66uJi8vLyprUygUkaNfbKsHDRJNfsPWmKaGw9G132+PW2+9lW9/+9uUlpZy8OBBRowYgcfjYe/evQwaNIgVK1YAMHPmTCZMmMDnn3/Oeeed1+NjUSgU8UG/EEcZaT5V4Pz+AC5X7/IhMzIyeOWVVwB48cUXAfjJT37CM888w8MPP8yrr77a/Nx169axaNEiEhISevWeCoUi9vQLcYwGdXV1zVbitddeyze+8Y1TyhPz8vKYN29erJaoUCjCiBLHLpKSktLKShw9evQpzxk/fnw0l6RQKCKICsgoFApFCJQ4KhQKRQiUOPYQIUTIxhHt3a9QKPoWShx7QMvGE20xTZN7770Xv98fg5UpFIpw0S8CMgcOQEPDqfc3Nmq4XNBZF7HERGgqfukSLRtPtMXpdHLhhRfyu9/9jh/+8Iddf1GFQhFX9AvLsaEBUlJOvSUlgcsl/w31uH0LJaw2L7/8MtOmTePBBx/kxhtv5O2332bZsmV4vV62bt3KmDFj+Ne//sWuXbu46KKL2LVrF3PmzOGNN96I3h9AoVCEnX5hOXaEroPfDw6HvHW3R+6NN97Igw8+yL333ktdXR0XXXQRJ06cwOl0MmPGDK677jpSUlLIzc3lqquuYsKECQAcO3YsAkejUCiiRb8XR5ACaZogBDid3RfIoU390FJSUjh+/Hir7t9f/vKXWbx4MXv37uWqq65qvt+yrLCsXaFQBCkqKmrVqiySxEQcPR7PWOAhoBgYBhwzDOOBSL6nroNlSSvS6ezcD9mS0tJSQFbJZGRkEAgEmh87++yz2b17Nzk5Odx6663N9ycnJ4dt7QpFPBFNgWr7vm3HI3i9XoSASIxjipXlmAW8ahjGGwAej+djj8ez3DCMrZF8U12X1qPfT5cCNTZpaWk8/vjjFBcX88ADD7B06VL27t3LqFGjALj88svJzs5ufn5JSYlqPqHol7QnUNEg1HgEj8fLkSNQXq53K6jaFWIijoZhbGlzlw7U9fT1EhOhLsRvnzwpBbDtNlqI4Bbb6ZS/3xHp6en86Ec/av75jDPO4JlnnuFXv/oVpmmSkJDANddc0/z4s88+y913393Tw1Eo4pZQAhUtcWw7HmHmzHz27ZMuM8sK/8C9mPscPR7P1cBqwzB2dvS8lm3QLctqlUfYXvvEhgYLXRft+hgtS8PhkKMU2ktL/Pvf/87+/ftZt25d80kwZMgQbrvtNvbu3ctPf/pTzj33XJKTk/H7/fj9fm6//XaGDBnSaa6jZVldau9eVVXV6XP6Cv3lWE7X45gyZUpzsxWXy8WUKVMiMqIgFGPGjOHll19m06ZNTJkyG7f7DKqrj6LrUFNzgvLy9ptQ94SYiqPH47kAuAD4r86em5OT0/z/mpqaVt1w2sPv9+NwuDoMwFiWtCJdrtCBmhtvvJEbb7zxlPvtJhN2OzMbl8vV5QYUuq63Oq6O6Orz+gL95VhOx+NYtGgRq1ationPEeDiixdx9tmLqKmRaXjBWIJGTs7AsL5XzMTR4/FcDswD7gQGezyekYZhbIr2Ouw/rs/XPT+kQnG6Eo35LaFobIRDh+Q2Oi0t8u8XEynweDwzgL8BXuBd4A1gQizWAkFB9PmkUCoUiviipgb275f/j1YiSKwCMluB1Fi8d3vYW2rbguzKeAWFQhFZLAuOHoXjx6UoRvO6jHlAJp7QNHnz+6Uf0uGAtWvf5Zvf/CZXXXVVc9DlnnvuCZnHuGzZMvLz88nIyIjB6hWK/oXfD4cPy/Le1NTuF2/0FuVhC4GuQyAgb+effwEjR47kjjvu4P7772fu3Ln8+Mc/Dvl7y5Yt6zdRTIUiltTXw9698hqMhTCCshzbpWXJYUsuu+wy7r77bh555BH8fj/79+/nRz/6EU6nk+3bt/P73/+eRYsWYZomK1aswOl0Mnv27FZ5kAqFIjRCyC10RYVsGBOi8VXU6Pfi+N57RaxfX8jcufnMnNm9CJsdybaTxm0yMzMZO3Ys1113Hdu3b+dPf/oTjz32GFOnTuW73/0uo0aNYvPmzTzyyCMkJCRw2WWXKXFUnJZ0p9TQNOHIEaitjZ212JJ+LY7vvbeZq6++sjmjftmylT0SSAhGsnVdJs66XC5+/vOfY5omx48fP+X3UlJSeOCBB8jMzKSysjIch6NQ9Cm6U2rY2AilpbBtWxHbthXi9eYzfXr004Va0q/FccOGda1KndavL+y2ONpomhTINWtWc/7553P33XfzySefsGfPHh5++GEAHA4HQgh27tzJfffdx3333cfUqVObR7oqFKcTXS01rKmBsjL46KMibr89KKZLlqyMqUD2a3GcM2deq1rMuXPzu/0a69cXcuDAfp599hlSUlKorz/JQw89gqZpfP/73ycjI4MPP/yQkpISZs+ezWOPPcaZZ57J1Vdfzf3338+cOXMoLS1lzZo1LFiwIAJHqVDEJ21rofPzW19/liV9i1VVMk1n61YpppZl4vf7KCoqVOIYKc47bybLlq3ssc8RYO7cfHbs2NXqPsuCxx9/qrk35EMPPQTA2LFjueWWW5qfd9NNNwGocQmK0xKv18vKlStD+hztNJ3GxqB/0euVYur3+3C53Hi93Tdmwkm/FkeA887z9ngr3R4tI9nt1WQrFIrQpYZ1dbIM0OmU9dE206d7WbJkJUVFyufYp1E12QpF9xACKivlVjo5OXSazvTp3m6LomVBTU34L8A+eUlrmhYXs6F7U5MthGg1bkGh6M8EAnIbffSobBoRrvxFnw/WrdN4773Ou3R1lz5pOSYkJFBVVUVGRkbMBaYnNdlCCKqqqkhISIjs4hSKOODkSbmNFiK83XRqa2HtWo1jxzT8/n7Y7LYn5OXlUVZWRmVlZYcW5MmTFg6HHrVidSGkNdnZFlvTNBISEshrr0uvQtEPEEJGosvLZbf9LrRg7TIVFbB2rY6uw6BBUnzDTZ8UR6fTybBhwzp93h/+cILjx9O5+GKLLjy911iWdDanpckPTPkhFacrpilFsWVT2nCxZw9s2KCTni59l5EYrgV91OfYHRITYc0anU8/jfz2W9elMFZXy2+yFkMKFYrThoYG2LcvaCiESxgtC95/X2PdOp2BAyPf17Hfi2NCAuTkwKZNGtu3a6c0kogEqany2+zAAemLVChOB4SQhsH+/VIQwylePh+sX6+xY4dGXh643eF77fbo9+II0teRlwfbtmm8954WlW7fSUny3337ZPslhaI/Y5qyBLCsTIpiOMWrrg7eflujtFRj8ODoNbztkz7HnuBwwJAhsHu3RmMjzJ4tIt4OKSFBvu+BA9J6zcyM7PspFLGgoUG6kSyr59Ho4uKikMnfR4/Cu+9qaJpGtOeZnTbiCNLUHzwYDh7UePddyM8XRDqbxq4CKC+XWwMVqFH0F+xtdHm5tBQ7m//eHsXFRdx006kNJ/btg3XrdNLSpKsq2px2l6mmQW4uHDum8c47WlS2vC0DNQcPtj8jW6HoK5imTOo+cqT32+iiotYNJzZtKmTHDo2CAhl4iYUwQj8Wx8OHoaSkfcN40CCoq9N46y2N2trorCk1VUaw9+2TibEKRV/k5Ek5wuDkyfBEo+2GEw6HA5fLTVrafLZti17gpT365bba74ebbnKxaZObRYssrrrKCtkcIisLqqs1Vq+GhQsF0ZiLlZgo17d/vwwSpadH/j0VinDQcoRBOJO67YYT69YV4nbPx+mcRU5O7Bu69EvLUdNg9mwLy9JYvtzB4sWOdrfP6emg6xqrVukcPRqd9blc0g9ZVia3JWpWtiLe8fmkS6iiQu6AwlntAjB6tJfhw+9m4MBZ5ObGXhihn4qj0wkPPmhy6601JCYK3n9f55e/dFJaGvr5AwZIv8lbb+mUlUVnjW39kCofUhGvnDghXUF+vzxnwy1cpaWwapWOpmkMHBje1+4N/VIcbSZP9nHnnQGGDROUl2s8/LATwwj9yaakSCvy7bd1DhyI3hpTU6Vz264oUCjiBTvocuiQ3Eb3NBrdHkLArl0a77yjM2CANFLiiZiJo8fjyfN4PH/yeDxbIvk+2dnwk58EmDnTorFR449/dPKPf+ght7KJidIPuXatzp49kVxVaxIS5O3gQZnXFQfd2BSnOfbcaLsEMNyJ16YJhqGxebPMXwy38IaDWAZk5gJvAGdH4sWLiop4551/M3r0JaSlweDBBXg857N162xWrXKwa1c1Z531PFOmzGTs2GDSaUKCFNR163R8PsGECdFRKqdTWpHHjsmk2tzc8Pt1FIrOsCx5DlZWSsHasSN0cnZvaGiADRs0Dh+WFS/xmvcbM3E0DOM1j8dzfiRe2x4J2dDgQ9cfQ9M0AgE/YAEXAH9jz55B7NlzFStXfpF587zMmiXnvezaVcCECfMZMcLL5s0afj9MmiSi4iDWNPkt3dAgv7WHDGndRl6hiCQNDXIb7ffLL+pt20InZ/eGmhpZ8XLypBTGeKbPpPKUl5d3+bkrVqzA5/MhhIlp2vtn2wJ8FzgXeB04B8taR0HBV1i//kI0TcM0A+i6g3PPvZEZM26gsHAmFRV+Jk/2RzWCFgjABx/oZGZaaFpV9N44wlRV9Y9j6U/HYVlQVaVx/LgDt9vC5ZJCuWbNqlbJ2WvWrGLEiHE9fq+KCp2NGxNwOiEtTVBdHZ5j8PmgoaGB8vLwOu37jDjmdKOw8rLLLuN3v/tdk+XoaBI9P0JYaJqOrh8GFmKazwDXA29gmj8FHgUEpmlSVPQXtm79G9df/wTbth0jPz+fr3xlZlS3AEJI38+JExqjRw+MeKljtOjOZxnP9IfjaGiAhoYcdB1GjGgdiV6w4FKef/6p5mmACxZcSnZ2drffQwj47DMNw9DIzQ02ZQkXjY1QW6uRkxPebVafEcfuYI+EfPTRoM9x164CUlMHUlt7jAkT5gOwYcNLbNjwEZZ1P/AImjYVIW4HGgBBINDIq69+H8sSrF3roKDgFr71ra/g8URnKpqmyW11dbVg3z7ZvCI9PT5ywBR9G8uSfsXSUieDB4cu0QvHNEC7B+OHH8rAS7j96ELAxo0a1dXhL6WJmTh6PJ75wFeBwR6P5x7gCcMwwlZU5/V6WbhwEn5/OgMH0iroYjN2rJc5c4p4992/U1x8DX7/jaSleamvvwAhSgENy7IQwsI0Tdas+RPr1i3h3nsfp6bmWNTGR7rdMg/zyBEZPYzESaY4fTh5UhYgBAKQnGx1WKLXk2mANo2NUrjsVmPh3nUdPw4vvujgww91XK4UDh/2hdWPGcuATAFQEMn3mDDBzwcfyBzH7OzQH87YsV7GjoWDBy1+9zudysoxDBiwixkzXmLYsACvvnoXgUBD06wagd/fyC9+cScgcLvd3Hvv4xw/HnmhtJPG7brW3NzIJOQq+i+mKVPFqqpkVkZKSvdq/NtrKxaKmho5/Kq+PvyBFyFg82aNV15xUF+vkZwsuOiiOvLywmsx9Mtttc2gQRZXXCH46CPZxzErq31/x7Bh8LOfBVi82EFJSSIbNnyd2283ueuuyWzc+BIbN76IZQUADSGkNenzNXLffXdiWQKn08EXv3gL11zzlYiKZFJSMDm3tlY20FBWpKIjhJDnypEj8ufU1O5/qbbXViwUZWUyV9jtludnW0pKipqzQkLt6DrixAn4618dbN0qLZ0pUyxuuMGkvt6Hpilx7BYJCeD1CoYPF2zcqHHihMagQaFPjgED4K67TF56CTZt0vnf/3VyzTWzuekmL7Nnf7XZbymtSR/2ttuyLHw+k1de+RNLly4JS8pDRzgcyopUdA2fT/ZbrKuTX6x2g2fbCpw4cRoLFlzS6eu0bStWVFR4yjkuBHz6qUzszswMbYiUlBTxxBOXEAj4cDrd3HXX6i4L5PbtGn/5i4OaGo2EBMENN5jMnSvw+SLTbb/fi6PN0KFw5ZWCrVuhpERus0NFf10uuPVWk8GDBUuXOli61MHhwxpf/aq3+UMcOnQyu3YVoOsDWbbsLoSQ224hBH6/j6VLl4Q9cTYULa3Imhrpi4xliydF/CDTc+Q2WqbOBB9raQW6XC7++tdVnZ6ndlsxO3Lt9ea3etw0YetWjZ07Ow687NpVQCBgp9n52LWroFNxbGiAv/9dp7BQlumMH29x661mSKs0nJw24ggy43/2bMGwYYJNm+Q866ysU5+nabBokUVuruC55xxs2qRTUQHf+Y5JWprtp5Qf6JAhk3nvvZcoLpbbbofDwWuvvUQgEIjKVtu2IhsbpRWZnQ0ZGfFbdaCILHb6V1mZFKyUlFN3FC2twECAkFZgW/9iR5Hr+npZ8WIYm6moKGDixPa3yxMmzMfpdGOaPhwOd3PmSHuUlGg895yDigoNp1PwH/9hcdFFVlTO79NKHEGeKCNHQna2RVGRxqFDWrt+u+nTBdnZAX73Oyeffabzy19qfO97AYYODT5n6lQvY8Z4mTbtq+j6WmpqDvDqq89jWWZUt9oJCdJqtB3ueXmRH12piC8aGuTnb2+h2/Ovt7QCnU7XKVZge/7FUJHrY8egoECjpGQzL7wgt8vLl5+6XW7pZ7zrrtWd+hwDAVi+XGf5ch0hNIYOFXz964GozJ+3Oe3E0SYlBS64QLB7tyyAT0kJ3RVkxAgZqPn97x3s26fz6187+cY3TCZPDtZcp6bC5Mleqqq8TJmykaVLl9DYGP2ttqbJtfj9cqhXWpq0JNVWu3/j98ucxepq+SXf2ZCrllbgxInTTjkfu+JfBNlJav16ndRUKCtrf7scys942WV3t7u+sjJ47jkHe/fqaJrgkktMrrrKinrg8bQVR5Bbz4kTBYMGCdav1zhyRAs5ACsjA/77v03+/GfYulXn6acd3HCDxYIFwQ7jSUlyi1tRMZuHH17Fli0v8dprL2Garbfa4apR7QiXi+YSsL17pesgMzN6Iy0V0cE0pSAeOxYsGOhqUM62Ao+G6PDcmX/RsmDHDo0PPpDXi9vd8Xa5q35GIaCwUOf//T8dn08jK0tw221m1Jq/tOW0FkebgQPhsssExcWyv1yoYE1CAnzjGyZvvin4178cvPqqg7IyuOEGqzkC6HbLoEh5+Syuv97L1VffxObNhRw61HKr3chTTz3EnXfeE/EE8sTEYGv748fl2sIx80MRWyxLprRUVMj/JyeH9zPtyL/Y0BBM7M7LC37hjh3rbXe73BU/Y00N/OUvDrZvlwfi9Vp8+ctmTF1DShybcLlg5kzB4MGCjRt1nM5T50zrOlx1lQzU/OUvDtaudXDkiMY3v2k2d89xuaS/b9s2jYkTZ/Gtb3l5//0ili5dgs/XiGVZbNiwhi1bNkQlgdy2KOyh68eOydyznuS6KWKLZcl8xYoKu7olcruBUP7F48elf7GhQWPIkFN/p2Wgsu39HfkZt2/XeOEFBydOyITum24yOffc2Dc1VeLYhhEjICvLau43l5Nz6gno9QoGDTL5wx8cfPKJziOPaHz3uwFyc+XjDodsN7Z7t0ZDA8yeLb+Jn3rqITZsWNOUFxlMII/GVtuOagcCMvXHTtBNTlYiGe/Yonj0qPz8EhNPDbZ0p3qlJ+zfL/2LSUmhE7s7I5RwNjbCa6/prF0rL7CJE2WKTqgMkligNlghSE2FCy8UTJkiKCsLXWI1dqzgZz8LMHSooKxM41e/cvLJJ0GV0XUpkKWlcj72xIle7rzzHtzuBBwOB7quNyWQB7faxcVFET82u6mupsnZHfZ4BtV9PP6wLLnd3LNHWv32Z+dsY9LY0eUnn7yfm25aFNbzyLLkLmjtWp3MzPCNMti7V+Ohh5ysXevA6RRcd53JD34QP8IIShzbxeGAadMEF15ocfKk3I62ZeBAOYJh2jSL+nqNp55ysHZt6z9pbi7U1Mj52OPGSQvyBz+4j/vvfwq3O6FZJDdsWBP2E7sjXK7WIrlnj7wQ1STE2GOacgv7+edSFO3Pqq0o2oSKLoeD+nrZmPbDD8M3Q9qyZIrOI484KCvTGDxY8NOfBrj44ujkLnaHTrfVmqZ9QQixLBqLiUeGDIHLLgvmRObmtt5mJybK5PClSwWrVzv4618dlJa2DtQMHCjnY69cqbFgQdCXM2HCpFZb7WhW19jYke1AQF6IdmJ8Wlr7F6MiMvh8MvpcVSUteTsDojM6iy73hKNHpX/RNEP7F3tCRQU8/7yDzz6TKrhwock113TcFSiWdOX0v0fTtLnAi0KIDyO9oHjEzon86CPZmy4jo3WCta7DF79oMXSo4MUXZaCmrEwGauw+eenp8pt49WqdefMsRo6UTu8777yHLVs24Pf7op7y0xJ7y2aa8iQ+elSuOT09dJmlIjwIId02x49L94aud98PHI6+iy3X8+mnGu+9pzFggExj6y1CwKZNsotOQ4NGerrg1ltNJk2Kb19OV8TxFmAf8DVN074NrBJCvBnZZcUfug5TpghycgSFhbIVU9umyLNmCXJyTBYvdrBzp86vfqVxxx3BiprkZClCa9fqTJ8umDRJtDqx20v56U1r+u7icEiRFEKmi9jtrbKy5JdEvG19+ip+vwyyHD8urXZ769xTetN30cbnkwURn30Wvsa0J07AkiUOiovliTNjhsVNN5m9OtZo0RVxdAAm0AjMBkZqmnYxsE4I8f8iubh4JDcXLr9cUFQkgy05Oa23n3agZvFiWVHz8MNObr/d5Jxz5Lek2w2DB0sL9MQJOO880XxiFxeHTvn5wQ9+gd/fGLWtNkjLxY6I+nxydrGuS0tywAApmCrK3T0sS1qJVVXSStS0yMyD7gnHj9P8pT9kSHg+2x07ZBed6mqNxETBl79sMmtWdIbVhYOuiOMSIA34F3C9EOJTAE3THgBOO3EEaQGef77gk09kJ5LMzNbb7KwsWVHzl7/Ali06ixc7ufJKkyuukE5nh0MK5N69GjU1MG+eICUluD1qm/Lz2GM/i1rKTyjcbnmzrcnqankMGRnSmlTb7vaxLJmyUl0t/3ZC9N5KDDd79sgWfUlJslCgtzQ2yi46BQXBLjq33WaestPqKr3p/9gbuiKOu4CvCyFO2HdomuYGMtv/lf6PrsuRrdnZobfZCQnwn/9pMmKEYOlSnTffdHDggMZtt5kkJclv5txcqKzUWLlSim129ql+SDkcTDbX7ajONRq0tCZNU0bw7ZZYth9WWZRSEBsa5La5uloKotMZfzmlfj9s3ZrAkSN6cxlgbykp0Xj+eQfl5RoOh+yi05tIdG/6P/aWrojjjUIIs+UdQggf8L3ILKlv0dE2W9Pg0ksthg0TPPusg/ffl37I73wn0Nw6PitLTk5bvVpj9myL0aNbO9gzMwfywAN3EQj4cTgcHDp0gOLiopgJpI3DEZypbZqy8YGd7jRggLSMEhJOn4i33y+3zLW1wbxRp1N+mcSjn/b4cVi3TuPIEQejRvV+jX4/vPmmzqpVwS46t98eYPjw3r1uT/o/hotOT922wqg4FXubvXOn3Ga3jWZPnmz7IZ2UlsqE8dtuC/ohU1Plt/a6dTrHjgnOOUe0crDn5g5jzZo3ee21l3j11eej0gKtOzgcweO1LCkOVVXygrO3kMnJ/as7kN8vj7OuTgpiICDvd7niz0JsiT0mdcsWjaQkGDiw9/mFBw7A8887OXhQQ9MEl15q8oUvhKeLTnf7P4aT0+R7PfLoOpx1Vvvb7JwcmTD+wgty/sXixU4WLTL5j/+QJ6cdqNm1S+PYMZg7VzRbZlOnzmDnzg8IBALNib7RzofsKrreOsgQCMitZWWl/Lm62tGcw+d2SzGJ925BpinFsLFRpmPV1UFFhZPa2mCeaDwEVTqjoQHee09jzx6tuWt8dXXPX880YdUqnTff1DFNjUGDZBedcePCl6LTWV02yHMsEl9GShzDTE6O3GZv3gwHDrTeZicmwje/afLvfxTgyy8AACAASURBVAv+8Q+dlSsd7N2r8Z//aTZ3y8nLk37IFSs08vOt5nrtlom+0e423hucztZb67o6QX29vCg1Lbj9tAXV7Q7+jsMRPQtMCHmxBwLy1tgobydPSmvYLq+0hTAlxeq0b2I8ceSI3EYHAhpDh/b+73rokEzo3rdPmp3nn2/yxS9aEQnOtdfQwu+XrhyHQzBpkj/s79vvxTEWNcNJSZCfL9i1S+aN2T44kCflxRdbjBwp+OMfZeOKBx6QCeP2N25WlrRQ3npL5kPm5NBBPmT0uo2HA4dD+iJbXkR2RLe+Pli+aF+8DocUTPvflqKp6/Jf+7mhLnghgjdb5GwR9Pnk/+1/Wz7f4ZDvlZAQfp9hpJtEtMQ04bXXNvP224VMmTKfSZN6936mKc/LZct0AgHZc/GWW0zOOit6F5odDNQ02a1/zBjB0aPh9/71a3HMzJSTyU6ebL9lfKTQdTjzTNlIt7BQo6JCbrPtC3jCBMG99wb4v/+T5VSPP+7g6qtlZE/Tgj664mKN1NRELrmEU/Ih23Ybj2UkuzfoevvpQJYV3Nba4tn2C8+2QNuKo32f/Xz7Z1tY7ZvbHb2gSXdGnPaWqipYsmQzTz99Kabpo7Cwd9He0lL485+D1uK8eRbXXWdG7dqyLOmeMU3pwpo4UZCYGLl+AHEYRwsfiYmCUaPkhVdbGxsrMjtbNtIdMkRw6JC8yG0yM+Uo2IsvNjFNjddec/CHPzioq5OPO52ytruyUmP5cp2DB+X9thX55S/f3tzlx+Vyk5k5kMWLH41a84poYAd1EhLkF1xKirTCW97s+1JSWt9aPpaaKuvF7fuSkuRrulzRjSZHqklESyxL+q6XL9f55JNCTLN1tLe7BAIyEv3gg0727dPJyhL8138FuPnm6AljVZV0DYwcKfjCFyzOPltE3M8bM8vR4/FcCFwDlAPCMIz7I/E+Tqccy3r8uKwZbjm7N1okJspE78GDpUM8OTnY+snphOuusxg/XvDnPzv44AO5zf76103Gj5dqnpkpSEiANWt0zjpLcPbZwWj2Ndfc1Jzy8+CDP2q2SKLRSFfRfSLRJKIltbWweXNwcNy0afN5++2eR3v37JFVLqWl0izPz5e+xWiJop0rOmyY4IILRFRbmsVEHD0eTzLwDDDJMIxGj8fzD4/Hs9AwjHci8X6aJv14iYnSkWw3DI0mmgbjxwsGDgw9r+bss4Pb7D17dB57zMEXvmBx2WVyz5CYKKPZu3drHDokR8zaSePTp3tZvPjRZosk2o10FV0nnE0iWmJZ8PnnsmGE201zJ52uRHttWlaiDBvm5Z//1HnnHZm3OGiQ4OabTSZOjM72q7FR+hUzMgQXXSTIy4t+elSsLMdZwD7DMBqbft4AXA5ERBxtkpPlWNYjR+Q3UncGEoWLrCxYtCj0vJrsbPjxj02WLROsXOngjTccfPKJxnXX6aSnSyHNzYUTJ2T7s7PPFpx1lsDhaG2RaJrW1Eg39lU1ilMJR5OIltjWYmmp1qrSpaXYdTTtz36uXYmi61eQmvoq1dUuNE1w8cUybzEaZaKBgBRFp1MwZw6MGiU6dXsUFxfx9turMM1L8XrD93eNlTjmACda/FzTdF+7lJeXd/tNqqqqQt7vdIJp6uzbp5OUZMUkz27sWEhMdLBlixuHA9LTg9/ICxfC8OEuXnklld27dZ54Ip0vfrGOs8/2NT8nMRE2btTZsUNw7rmNjBgxjt///lWKizeRnp7Jk0/eRyDgR9cdfP75p6xZs5qpU2dE/0DbUFNTE+slhIVIHsf27VspLt7E9OmzOvzMhIA9e5y8/74bl0uQni44eVIGIPfu3cIf/3g1gYAfp9PFN7/5OqNGnXvKa5w8Wc/evVt4661H8fszgScxzRuoroYhQwJcf30tw4aZNDTIPMlIIQRUVcn+kRMn+hk3zo/bHcyPbY/t27dyxx034Pf7ePHF3/LKK68wY0Z4zvNYiWM5spmFzYCm+9olp4cV8e39Xl6ebARgd1qORfVGdrYUyY0b5Ta7ZU7kuefCxImyecUHH+gsWZLG7t0WN94YHOaVlSW7dxcVpTF1qmDevEtYsOASAGbMmMnSpUt47bWX+Oc/X2bFitfiZnud3dMOBHFGJI6juLiI7373S51Gs48fl9ZiRYXG8OGnnr+lpQaBgL8pECN/njZt4Smvs3fvFp555osEArcBrwPpQB0XXHCEG24YjsOREvZjbEtNjbR+zzxTMHWq6FZTjp07P8Dv9zftkvzs2LGDRYsWhWVdsYpWbwJGejwe21CfAyyP9iLS0uRALbuVVCyw59VMny4oL5cnScv13XGHyTXX1OJ2C957T+cXv3Dy0UdBX8CAAXKr/dFHGm++Kf2RILduQ4aMOKWqpr9Fs/sbnUWz/X45rW/5cr25vVioL3a77E7XHR0GYrZu3UcgUAA8DaSTmrqJb3/7I268cXiXdlQlJUWsWPFrSkq6f041NMj0oMREwWWXWcye3T1hhKA7SdcduN1u8vPDF+CKieVoGEa9x+P5NvC0x+OpALZHKhjTGQkJQT/kiROxGVlqd/jJzRWsW6dRXi59kXaS8+zZjUyfnsBzz8lgzW9/62yVY+ZwSIGsr9d4+22NUaOk2LZXVaOCNPFLe9FsIeTUyKIijd27N3PkSAFnnjmf9PTWn2FLP2NHgZgTJ+D11x1s2vTNpnv243D8iO9+93uMHevp0lp72jHHrmxJTBScf75g2LDQ6VRdSZafPt3Liy+u5O23V3Hjjf3D54hhGP8G/h2r92+J3V/R7ZYfWiTnAXdEdrYsPdy2TQZrBg4MRtVzc+Huu03eekuwbJnOunU6H32k8dWvmkyeLP2VyckyVamsTGPZMo1p02bxwgsrMYz2u4wrgYwtbQUgVDS7tlYWA+zdq1FVVcSf/ywFacWK1oIUSqzaBmICASgokBUu9fWyrdi55x5k0KBXmDTpe91KEO9uxxzLkteXELKyZfx40W5aXXeS5adP9zJ48HjOPXdgl9feFfp1hUx30DSaI8ex9EO63TBzpkwaLyrSm5sbgBTsRYsspk61misVnnrKyaxZFtdfbzZbvQMHym/nDz7QSEiYxRVXeCkvD91lXFmQsaM9AbBvfj98+KHG9u0aLpdMz/ngg/YFqTOx+vBDjb/+tZGjR6Uf8ayzLK64oprx4/OAH3Z7/e11zAnVnLaqSrquxo+XI4+Tkzu2DEO5F6J9nipxbENamhSo0lJZrtay9Vg0GT5ctpPaskXj44/15i42IJPa/+d/ZAOLZct0Nm3S2bFD4/rrTbxe2Ybe5ZLWZkODRkGBRlbWLJ56aiVLlsR22qEiSHsCYFlw8KAsGPD55A7CtrA6auHV3mOHD8Pf/+5gxw4decmX4HDczZVX3smgQWf1eP2hcijbWq/f+c5qsrK8DBnSOom7M8sw0snyXUGJYwgSEmSgpqwsdvmQIIU5P1+QluZj1y55gdgnl8MhG+mefbbFkiUOdu3Sef55Jxs3yoi23Uw3MVFaHLW1GpWVs5k1617ee28DgYDyQ8aaUAJQXi57gh49qpGVJUtMW9JRUnfbx7Kzvbz0ks769TqWpeF0NhII3Af8FiEC7N59dq/E0X7Plmtoab0GAj4++6yA//mf8xg8uPU11JllGKlk+e6gxLEd7LLDo0dlrlWsJu9pGowYEeCMMyzee0/j4MHWib55ebI+e9Mmi7//XU49vP9+jYsusrj8cqvZZ2nXF1dXz+Kmm1Zz/HgBsJ9//lP5IWNFSwE466x8ampmsWqVRloaHc6Kbq+Fl/3YkCFe/v1vnbfe0mls1NB1wfz5JlOnbuOZZ57GNAMRaxw7bpy0XgMBaRHecsu8kMfSFcsw3Mny3UWJYwdoGgwaJC3Jw4elFRaO7sY9ITVVdhv//HO53XK5glakjGgLpkwJ8PrrDtav11i1ykFRkc6115qcd16wyiA9HWbO9FJT4+WTT4pwOJYghPJDxopx47w0NMyipER25u7N5D+fTwZbVqzQqa2VLzJtmsW119o7iRmnWJ3V1eFJZrdHZaSne3nooVWUlRUwZ07HUeZwWYaR6sqjxLELDBggRbG0VH4QsZq2p+swbpwgL0+EtCLT0uDmm03mztV4+WVZAfTcc07WrLG4/nqruV+kpkmR9Hq9JCWtZuXKB/j88zWq3DCK1NTAxx/LGdFut9wB9HRn4vfDhg06y5frVFVJURw3zuKaa6zm5iU2HVmdPUEImZDe2AhnnCFnsaekzARmdvq7vbUMhZBBHiEgKyv8CqnEsYskJcl8yEOHZJv8lMgXDrRLaipccIFg716ZGA4yQm1bHGPGCH76U7nVXrpU5kb++tc606dbXH21SV5e8LWmTfOSmvpznnhiQ1NdrZu8vPkx6YHZ3ykuLuLddwvJzp6PwzELt1sGzXoriitX6lRWyg9/+HDBVVeZTJ0a+fnQ1dXyWhgxQjBtmiAjI7Lv1xKfTwpyRoY89ysrw98QQ4ljN3C5ZBT5yBH5zR+LhHEbTYPRoyE312LrVjkXpOX8bF2HOXMEM2YEWL1a+p+Ki3Xef19jzhzBFVeYzdvylo78MWPms3s3rF79GLNm5XP55TNbNelVdB8hYM2aIr7znUXNUdwf/GA148f3zGpqbJTD2FavDlqKQ4YIrrzSZPr0zhs19JYTJ+T5n5cnyM8XPZ5H3RMsS2aR2MUbkeyupcSxm9hzXhISZH/IWCWM2yQny16RY8YINm+WNdrZ2cE1JSbCVVdZ5OdbLFvmYMMGjXXrdDZt0pg/3+LSSy0yMoLbrZKSYJLxmjUOVq78Gl7vVxg1CvbvL2DePJXy01UCAemK2bFD44031rXKQfz004Jui2NtLbz7rs6aNUGf4tChUhTPOSfyolhXJ63F7GzBxRcLcnOj+6VZXy+/aHJzpasr0u+txLEH2P0h3W65zXa7Yz92dOhQuPJKwYcfyjrrpCRabXMyM+GWW0wuuQT++U85AfGddxwUFurk58vxDFlZbROJTbZseZbi4hcADcsKsHixm9/8ZiULF3pj5nuNd6qrYd8+jZ07NXw+eSGfe+58Cgt71nT2yBF45x2dDRt0fD6pCKNHW0yf/gmm+ToZGfPR9ch9YdXXyyTujAzBggWyaXM0Mzf8flmHPWCADJBGq1m1EsdekJoqTftYJ4zbuFxwzjmCUaMEW7bI/n4tSxBBWr3f+pbJwYMmb77poLhYiuTatTqzZgkmTrwUp/NXBAJyPg0ILMvfNGJCEAj4+Oc/11FVNZvBgwVjxwoGDVL+ycZGmRe7a5fslKPr8gvJzm7orOls26oSy4I1a3ZSUJBAWdl4QIrimDGVXH31AJzOTTz5pLTwly/v3WyY9rAtxfR0WQM9dGh0RVEIuQanU7qzon19KXHsJfGSMN6SzEzZ6efAASmS1dWtqywAhg2Db39biuSKFQ4MQ2P9ep316z2MHXuA5OQX+PjjnyFEAE1zNDXQlflxHs98Bg2S/ffWrZMHm50tGDMGBg0SpKfH/m8QDQIBadXt2aOxf7+GEPLzb2+72V6kuGVVicORy/z5hWzbNpjKyilNz2hA015BiKc4cGA3Ltfqbtc1d4faWulXjJUogvyy8ftlSW9GRmxyjJU4hoF4SRhvia5Lq3bwYMHOndLvZedGtrxwhw2Db3zD5Kqr5MjNTZt0SkoGAncxaNC3GDq0gAUL0nG7LXbtKiA1dWCrIU27dhVwxhnzSU31smULCKGRkCAYPlwGCbKyYhvZDzcnT8rP+YMPEjhxQsc0pdXc0s/bXXbuLMDv9wJfIxC4jnfesc3wA8Afgf9DiKOAwDQdzRZme2WEPUEIGWSpq5NfdAsWiF6lF/UU05S7sJQUmfMZS9eNEscwEU8J4y1xu2HqVMHo0YL335edXVJSZJ5jS3Jz4atftbjqKou1a3UKC3UqKlKoqLiMnTsFM2dajB6dwpIlczFNH7ourUnTDJzSrsrv1zhwAD77TKpwYqJsS5WXJwgEdLKyYv/l0VX8fli3rojCwnVkZ+eTmTmrafa1zuDBUhBLSorYuLFr2+WWlJdDUZFOYeFdwD3N948efZyzzz7Mm2/OxLJONlnu7mbL3X6trs6G6QghZKecxkb5ZTZnjpyTHgvL3w64DB4sc3ZjvftQ4hhmWiaMm2b0B3m1R1qajGpPnCgoLpbT6dLTT7XqBgygebCXYWgUFOh89plOQYGDgoLpyD7FL2CarwJHkNZM622dy9W6Jtjnk2JZUqJx4kQiaWk6mZnSMsnKEs2ljbEO8FiW3FLW1MCxY7Jx8AcfbOb55xdhmjIF54YbnqCu7hhDh3oYNmxhhz0NQz02cKCXrVt1tmzRKCmxvyGSSEtrZMiQLVxwQQIzZpwNjGfChBXN4gecIoS9Sej2++Uu58QJnalTBWecIRgY3o5f3VpLY6M8R6MZcOmMOFlG/8JOGD98WF5s3e1uHEkGDYKLLxaUlgZF0p7M2BKnE7xegddrUlpqsm6dzsaNFidPTgN+AzyOpq0FXgWWU1l5gJKSopAXa8tofmKiRVqa7Bb06acQCGit3jMzU/osMzMhKUnOJk5MlL/vcvXemggE5IXY0EDTvBWNY8ekUFRXa+zdW8SePQWMGzefM8/0UlFR0Dz3ORBo5JVXvo9lCZxOF3fd9VaHvr/gY8MIBK7h2WdHUlnpRAh5EAkJgnPOEXi9FmeeqaPrratK2opfOHyK9fVS/J1OmDxZkJFxkhEjYhNJbBlwGTYs9gHNtihxjBAul/zAKypkGkQ8+CFtNE2ubcgQwcGDgq1bNY4fl0nkoSzdoUPhS1+yuPZaeOutj9m82cGRI+OxrIXAQoQIUFBQwLp1r/ONbyQwY8Y5Hb6/rssLoe3FEAhIsaquhs8+k/7LtmKYkCBneCckyL+xyyUvLk2TryuEtACFkFUUthhKQdRa1eEKIX8vMVG+XnV1ES+8IC29ggJp6bX07cl0JgshLAIB2vX9maYM0pSV3YwQXwCmNm9fnU7B5MkWHo8cTB8Na9k0ZdS5sVEGWWbPlm4OlwuOHo3OqNW2NDbKz2fgQOLWzaLEMYLYo1QTE2U0O178kDa6LiPtw4YJ9u8XbNumUVmpkZER+lvc5YLLLx/P5ZdDXZ1JcbHFqlX7KC8fASzEshbyzDNSdCdPtpg0Sab6tBSAvXu3UFpqhPSTOZ3y1l4ARwhpaVqWtPzq6oJNB4SgKd1ICp5903XpF0xI6Dxh/9NPT7UCL7vs7mbfXmrqQF599a4mIXQ1H8MPf7gaw/iQEycm8eKL2Rw7ptHY6ARGACNwOhsZN66GuXMzmxu9RgPbStQ0WVI6bpxoVWYaC6JZ4dJblDhGgfT0YAPdQCD+cgJ1HUaNkjWypaVSJA8d0khNlT7IUKSkSB/mkCGHePzxfEzzIuALuFxf4NAhJ4cOOXjrLXA4LAYPPkFKynaGDq2msPDrmObRbs0csbGb+EaK9iLALbe3Q4dO5uOP1+N2n8/nn5/HW29pfPbZXGpq5rV6rczMes45J4Fp0wTjx+u4XBlA5K00n09aiaYprUSvV1bRxIMInTwpxTEnhz6R7qXEMUokJUkBiqd8yLbouky2HTpUUFYm+PBDjcOHZdeYjIzQVtfYsV5+9KNXm7aYgzDNIlav3syOHRZCLMA0z+HgwXRgHrt2ARwCduP3GyxffpKLL9YYNqz7U+ciQagIcF0dlJbKhPoDBzT275/DwYNzMc3WH15CQi2NjcuBd9C0NZx//tdOmd8SKQIBaSE2Nkq3w8SJMHKkbAQRD+eYdJcEAy7xtHvqCCWOUcTplLlbdgAg1nXZ7aHrcp1DhgiOHRN89plGSYm8yjIyTo0qt6zL/u1vL8bvbyBoJWUC+cjpu7MBDzARmMiOHbBjh3zWgAGCjIxqdL2E0aMzmTBhBFlZMpqdlhY5n5RlSWE5flzj6FE4enQ2R4/O4aOP4MgRmUDfFk0T5OYGGDdOZ+xY2Qqurm4bTz55a7PVmZo6kBUrft2rNJuO8PtlonZjozyvRo+WlVGDBsWP/04IuYXWNHk+paV1/jvxhBLHKKPr8tszMVFGs12u2KewdMTAgTBwoByKtH+/xscfQ2WlRmKi3HK3FHc7OhsURg2How5NW4Vl/QuHw83llz9CVVUKbvd8GhpGc+CAtMpqajRqajKAGezdC+++G3xdh0MK5IABkJoqmoI5ojmQYgdl7FG2IEWvvHwfR47sIyNjLMnJQ2lo0KivlxfsiRMaVVV+6uqcdDS+3e0WDB4sq0RGjBAMGyYYOVLQ2FhDenpLn4P3FN9kd0eWdkZjoxRy05THPGqUYMQIKYjx9iVrtxTLzJTnULytrysocYwR9iCvw4dj3x+yKyQnw8SJgjPOgIoKQUmJTCi3S+ZSU1v77DTNwdy5X2PWrJsAWsw1OauFqMhoimXB0qXPsHr1KmAMMI68vEW4XOM5flzOv6mqklF/u8a4a4xpurWH3S3kCLm5KQwZkszAgbIFV16eICdHBjBCWWKNjafeZ1vQK1b8OiylfaYpXTANDTT9nQVnnSWrntpbV6yxLHk+JybKYF+8+de7gxLHGJKQIH18FRXSiR5P6T7tYUfgc3MFM2YIyspg927ZKi0lxcu3vrWaAwcKmDjx1OFPQMi2/LoO55xzNmvW3I1prsbhcPO1r01n7NjRgNxC1tRATY1GXZ28+OrrNQ4f3kdBwd+wLCealsDgwVNITR3EZ59txLJ8gLxp2kkmTZrHOeecT1KSFPoPPvgz7777IFCGrlvMnv2LsPkIe1raZ5ry2Ozu1k6nFMIRI6Rgx4NftiNOnpTHYAdc4v1c7oyoi6PH49GB/wQeBBYYhvFhtNcQTzgcslNOcnJwXnY8b7NbYqdjjBwpqK8XlJdDSclMcnO9ze3zU1NDO+DbltV1VA7ncsmtWVXVJvbvDz6+YsXLCPELwEQIOHRIQ9cdCGFhW6Wg4XQmcsUVFzJ2bDBanJg4kfXrj2KaVrsC1lHpX0d0pbTPzsOsrw9aofJcEEyeLH2tmZl9Q2D6asClM2JhOU4DNgP1MXjvuGXAgGBddl/YZrclOVlG40eNEjQ2Co4ehYMHNfbvl+WDdrK1aZ5aVvelLz1Bbe0xJkyY32y9tRWmUKV4toXWsr2aEBa6rjclkAe39m0FqistxNorC+wKLdN/LCtYlePzBZPPU1MFI0dCTo6MLEcy8BQJ7AoXh0MWFfS1c7Yzoi6OhmFsA/B4PNF+67jHbn9mV9XEazS7MxISZFXN0KGC886DmhpBZaVGaancghtGy9nGjfz1r99HCNEsQsApwhSqTM9O0N648SU2bnyxuTFDS7FtWefcndrknrQEs2uEfT55s9F12Sh21CjZ0s3OH+3LFpbdUkxmFPQtUe8qERFHj8ezGsgN8dDPDcNYFon37C/YPj17m90XxbEl9qTD9HTB6NEwceJJJkyYx7p1cmYxaAhhl+T5MIwCHA5aieeyZQ8wY8Y1HSZoz5791S5ZgbrevjXZkrZ+w3Hj5tPQIAUhEIDKSr05UKJp8t/ERGkBZmbSNM9HCmFSUv8RD3sLHQ8txSJNRMTRMIxLwv2a5eXl3f6dKhne7LMkJUFFhU5FRW3zRdjXqampYcKEcfzhD69SXLyJ9PRMnnzyPvx+P06ni/PPP4f6eo01a1z4/XKb/PHHa9i9ewMXXvhLTp6sZMSIuWjaJPbuPdFcHpicPInp0yehaXDkyAk0Lehf3Lr1rVajHwoKnmXDhpe4+ebXGTLk3KYWZG1ruSfzpS+9zr59GxgxYg4JCZOoqTlBSoogM9OisnItO3duY/r0WZxzzgwSEkTIbjKyuUVU/rQ9oqama3OrZc6ihsMB2dkWLpegujrCi+sGkbjW+0y0OicnJ6q/Fy8MHQq7dwsgm4SE2M+qCQfZ2dksWHAJCxbI79AZM2aeMtx9+vRVPPXUQ2zYsKap2YOP0aNr+dKX7mH79s0YxmImTsxn1CgvPp8UN7/fFrnWtylTLmLt2ifajH7wUVe3halTL8DtprmZhcMhmhtauFzn43afj9vd2vIrLi7iN7+5Gb/fz2uvuVmyZGWfHjqW3cn4QLvsb/jw2HXl7grhvtZjEa3OBO4A0oFveDyelw3DKIr2OvoKmibbeA0YEGyBFo+lh70h1HD36dO93HnnPWzZsgG/34fL5WbhwnwqKor47/9ehM/nw+2WwjRzZmfCNJPZs1eydOkSXnvtJUwzgMvl5rrr5jF9evfrnYuKCvH7/ViWid/vo6iosE+LY3vYidwDBshsgf7wxdwdYhGQOQ481HRTdJHERJk2c+yYvCUl9W2HfleYPt3LkiUrW1mVixc/is/n67Yw2QJ8zTU3nWKldhevNx+Xy0UgAC6XG683v0evE6+YprQW7QBhX07k7g19ZlutCJYepqZKK9Lnk4Gb/mRFtqWtVen15uN2u5utye4KUygrtSdr+sMf/sbOnR/0SmTjDbudmD2bPR5GFcQSJY59ELvTeGWlvPUXX2RXCGVNxoKpU2c0+0z7OkIEq3Kys2V2QV/PkggHShz7KA5H0IqM5zZokSAc1p9C0tAgz52MDJmv2N9dNd1BiWMfx7Yiq6rkyNC+VH4YTxQXF8XcGo0mdrJ6QoKsbFLnzKkocewH6DrN86ErKqQlkJSktkZdpbi4iJtuah0B768CaYtiSoocgVpTYylhbIc4zVhS9AS7bG/wYHkR1NUF56oo2qeoqPCUCHh/w+eTzXGdTntuUHzPb4kHlOXYz9A0GWVMTpZdcSor5QWhLoT26W0EPJ6xLcXkZJnEHW/jT+MZJY79FFnmJRN4jx6VVsPpFNXuDvESAQ8ndmOIpCSZlnO65ir2BiWO/Ry3WzYIqK8P+iMTEwlZB3w60x8i4ELI6HMgIHcPgwerHUNvUJfIaUJysvQ11dZKkTx5su+2RFO0xp7jLYTcKYQagqboPkocTyNs/SNfugAAB7BJREFUf2RKitxmV1TIC0tFtvsmgYAURV2XLdL6eo/IeEOJ42mIrssqiLQ0JZJ9DSGkPzEQkC6TvDxZCBCvnXL6MkocT2PaiuTRo7LpgD3uVBE/2FYiyM8rI0P6E0+HiqhYocRR0Uok6+pk15/aWlVtE2taWolOpywXTUtTwbRoof7MimZ0XV58qakyYCNnRkvrpD+1+o937DQce8REWpqyEmOBEkfFKWiajGQnJ8sE4poaWbttWdLPpXIlw0/LoVwpKdJKTE5WX0ixRImjokPcbplMnpUlcyWPH5dbb5BbbrXF6xn23GrbQkxKkn/j5GT1N40X1Meg6BK6Lrfbqanygq6tldZkQwPNM6lVpLtj7PnVpin/ZraFqJLy4xP1kSi6jcsVHD9qN7ioqpJ+Sk2T1qaKdre2DkEKoO3TTUxUW+Z4R4mjolfYU/uysqQQnDwpfZS1tfJxp1MK5elgVdpiGAgE51mnpMjhVImJ8u+ggip9ByWOirBhB2vS0+XWsaFBWpW1tcHytoYGrTk1pS8jhLQI7XGwuh4Uw+RkKYZutxLDvkwfP0UV8YrDIYUiJQVycqSI+HxgWRYQ7DUphBRK+xZvYiKEFD/TbG0RNjRozWV7thDG4/oVPUeJoyIquFzylpUlyMmRwQnb8mpokDd7yBME/3U4pFVm/xtuP51lyZtpBv8VQr6PLYQJCTKabI/DdblgwACT3NzwrkURXyhxVMQEXQ/6K1NTg/fbFpp9swMadvmcbW22tNC60+28rWXncASbATudQQvQ6ZSPORyhrUFlIfZ/lDgq4gpbkNorWxRCWngt/215a4umBYVM04K+wUhYoYr+RdTF0ePx/AaoB2qBacB/GYZRFu11KPommnZ6RL4VsScW3511hmH8zDCMh4FtwM9isAaFQqHokKiLo2EY97R5/9por0GhUCg6IyLbao/HsxoIFcv7uWEYy5qekwFcDFzbldcsLy/v9jqqqqq6/TvxSH85Dug/x6KOI76IxHFERBwNw7iko8c9Hk86sBi4zTCMyq68Zk5OTo/W0tPfizf6y3FA/zkWdRzxRbiPI+rbao/Hkw38AfhvwzD2eDyeLlmOCoVCEU1ikcrzVtP7/tXj8QCcAP7R2S/dd999EV6WQqFQBNFEdzJoFQqF4jRBpcEqFApFCJQ4KhQKRQiUOCoUCkUIlDgqFApFCPpF4wmPx3MhcA1QDgjDMO5v83gi8DhQCowHHjEMY3fUF9oJXTiOu4E8oAyYgUyq3xn1hXZCZ8fR4nlfAZYAaYZhxF2lVBc+Dw34XtOPo4AMwzBui+oiu0AXjmM08vrYApwNvGwXa8QTHo8nD3gImGYYxrkhHteBXyGr7kYCzxmGUdTT9+vzlqPH40kGngF+YBjGL4CpHo9nYZun/Rewv6me+zfAc9FdZed08ThSgR8ahvFrZPrTY9FdZed08TjweDxnAmdFeXldpovHcRNQZRjG04Zh/BD4bZSX2SldPI4fA+sNw3gE+DXwRHRX2WXmAm8A7TWMux4YYBjGQ8DdwIsej6fHbUr6vDgCs4B9hmE0Nv28Abi8zXMuBzYBGIaxA5jm8XgGRG+JXaLT4zAM417DMOzcq3itS+/0OJou2B8DIS3KOKEr59VXgCyPx/N9j8djWyzxRleO4wgwqOn/g4CtUVpbtzAM4zVkXnR7tLzOK4EGYFJP368/iGMOrf9gNU33dfc5sabLa/R4PG7gFuCeUI/HmK4cxy+BBw3D8EVtVd2nK8cxEmmpPA28AKzqjaUSIbpyHE8CMz0ez5PAz4E/R2lt4Sas13l/EMdyIK3FzwOa7uvuc2JNl9bYJIz/C/zMMIySKK2tO3R4HB6PZziQCVzv8Xh+0nT3Dz1N5VJxRFc+jxpgM0CTD3sAMDwqq+s6XTmOF4A/NbkGrgb+5vF4sqKzvLAS1uu8P4jjJmCkx+Oxe0fPAZZ7PJ6sFlvn5cjtBR6PZwrwgWEYNdFfaod0ehwejycJ+CPwpGEYW+O0Lr3D4zAM44BhGF8zDOORJh8XyOMxYrPcdunKefUOMAag6T4HMlgWT3TlOIYDh5v+fxyw6CPa4PF4Ujwej+0SaHmdZwGJwEc9fe1+UT7o8XguAr4IVAB+wzDu93g8jwKVhmE80iQqjyNPgHHAr+I0Wt3ZcSwFJgOHmn4lJVTULtZ0dhxNzxkEfBN4sOn2R8MwSmO15lB04fNIBx4F9gFjgX8YhrEidisOTReOYy4yaFkMjAa2GobxTOxWHBqPxzMfuBm4FLl7egK4DZhiGMa3mqLVDyMnDYwAnu1NtLpfiKNCoVCEmz5hOisUCkW0UeKoUCgUIVDiqFAoFCFQ4qhQKBQhUOKoUCgUIVDiqFAoFCFQ4qhQKBQhUOKo6BdomvY/mqad1DRtvqZpP9Q0bYWmaWfEel2KvotKAlf0GzRNuxdZGeEDfiSEOBnjJSn6MMpyVPQnfgnMAz5WwqjoLcpyVPQbNE27FtmV5cfAFUKIz2O8JEUfRlmOin6Bpmm3AT8BViKbtb6uadr82K5K0ZdRlqNCoVCEQFmOCoVCEQIljgqFQhECJY4KhUIRAiWOCoVCEQIljgqFQhECJY4KhUIRAiWOCoVCEQIljgqFQhGC/w/EeWGm0cPCIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "train_and_test_approximate_gp(gpytorch.mlls.VariationalELBO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objective Funtion 2) The Predictive Log Likelihood\n",
    "\n",
    "The **predictive log likelihood** is an alternative to the variational ELBO that was proposed in [Jankowiak et al., 2019](https://arxiv.org/abs/1910.07123).\n",
    "It typically produces predictive variances than the `gpytorch.mlls.VariationalELBO` objective.\n",
    "\n",
    "\\begin{align}\n",
    "\\mathcal{L}_\\text{PLL} &= \\mathbb{E}_{p_\\text{data}( y, \\mathbf x )} \\left[ \\log p( y \\! \\mid \\! \\mathbf x)  \\right] - \\beta \\: \\text{KL} \\left[ q( \\mathbf u) \\Vert p( \\mathbf u) \\right] \\\\ \n",
    "  &\\approx \\sum_{i=1}^N \\log \\mathbb{E}_{q(\\mathbf u)} \\left[ \\int p( y_i \\! \\mid \\! f_i) p(f_i \\! \\mid \\! \\mathbf u) \\: d f_i \\right] - \\beta \\: \\text{KL} \\left[ q( \\mathbf u) \\Vert p( \\mathbf u) \\right] \n",
    "\\end{align}\n",
    "\n",
    "Note that this objective is *very similar* to the variational ELBO.\n",
    "The only difference is that the $\\log$ occurs *outside* the expectation $\\mathbb E_{q(\\mathbf u)}$.\n",
    "This difference results in very different predictive performance.\n",
    "\n",
    "### How to use the predictive log likelihood in GPyTorch\n",
    "\n",
    "In GPyTorch, this objective function is available as [gpytorch.mlls.PredictiveLogLikelihood](http://gpytorch.readthedocs.io/en/latest/marginal_log_likelihoods.html#predictiveloglikelihood)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAADQCAYAAACUXaMkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deXxU5b3/38+ZyWRfIBCWAAEi+xKVAw4CoYpa0SoK9dZaS2sXrbfV3lZv9f6uu9zaunWx11q1tte11p0KARSUIBLgAAKKBAiyhZ0QspFZznl+fzyZhGWSDCHJTJLn/XrNC3LmzDnPmZnzme/z3R4hpUSj0Wg0J2NEewAajUYTi2hx1Gg0mjBocdRoNJowaHHUaDSaMGhx1Gg0mjBocdRoNJowuKNxUtM0DeBfwErAA+QCP7As63g0xqPRaDSnEk3LcYVlWQ9ZlnUPkATMjOJYNBqN5iSiYjlaluUAcwBM03QD/YDiaIxFo9FowhEVcQxhmubXgV8A71uWZTW23913363LeDQaTZvwm9/8RoTbHlVxtCxrIbDQNM0XTdP8d8uynm5s3wcffPCMj3/w4EGysrLOZogxQWe5Dug816KvI7Zo6XXcf//9jT4XFZ+jaZojTdO88oRNXwGDozEWjUajCUe0LEcf8EPTNM8D4oARwO1RGotGo9GcRrQCMiWcRXQ6GAyyf/9+fD4fTXUVchyHioqKlp6mzRBCEB8fT+/evXG7o+rZ0Gg0jdAh78z9+/eTnJxMdnY2QoT1pQIQCASIi4trx5FFhpSS8vJy9u/fT79+/aI9HI1GE4YOWSHj8/nIyMhoUhhjGSEEGRkZ+Hy+aA9Fo9E0QocURyll1IVRSnnalN5xnIhfL4Ro0iWg0Wgi48gR+Oyz1p8hdkhxjDaBQIB7770X27bZvXs3s2fP5qc//SmPPPIIa9eujfbwNJouQ0kJPPecYNUqT6sfu0P6HE/lnXcEhw6dbknathuXq3n979lTcu21kVtxf/zjH7n88stxu90sW7aMiy++mNmzZwMwa9Ys3nnnncgHr9FozhgpYdUqWLjQID0djh1r/XN0CnE8dEiQnX369mBQEkkwuLRUAOHF8dixY/z7v/87Q4YMYcuWLQQCAXbt2sXPf/5zqqqq6oWwd+/eXH755WRmZvL5558zevTos7gijUbTGMEgLF4sWLFCEIpntoU46ml1M7zwwgtMmzaNBx54gBEjRnDVVVdRU1OD2+0mJSWFq666iquuuorLL78cgF69elFSUhLlUWs0nZPaWnjrLUFRkWDgQIiPb7tzaXFshpKSEgYMGABA//79AbBtu9H93W43fr+/Xcam0XQlqqvh9dcF27YJHKeIhQt/S0lJUZudr1NMq9uS3Nxcdu7cCcCePXsYMGAAycnJje5/7Ngxevfu3V7D02i6BBUV8NprgqNHBYFAEU8++XWCQT9ut4fbbltIQsKoVj9npxDHnj1lnd/wZGxb4HJF9vrGuOmmm7j11lspLS2tF0fTNNmxYwc9e/Zk/vz5AFxwwQUMGzaM7du3M2HChBZfi0ajOZnycnjlFUF1tYotzJ+/lGDQj5Q2tu1n69aljBmjxTEsKtJ8usAFAkHi4s4uHzIjI4PXXnsNgBdffBGAu+++m2eeeYZHHnmEf/zjH/X7Llu2jOnTpxPflo4QjaYLcfQovPqq4PhxQZ8+atuwYVNxuz3Yth+Xy8OQIVPb5NydQhzbg+rq6norcdasWdx8882nlSf27t2bKVOmRGuIGk2nImQx1tYKevVq2J6b6+WOOxZSXLyUYcOm0q+fl717K1v9/FocIyQ5OfkkK3HQoEGn7TNkyJD2HJJG02kJ+RhPFcYQublecnO9ALRVFa6OVms0mpiiuhr+8Q9BRUV4YWwvtDhqNJqYobYW/vlPQVmZoG/f6I5Fi2MLCdd4oqntGo2maQIBeO89wd694Sve2hstji3gxMYTp2LbNvfeey+BQCAKI9NoOia2DQUFguJiQV2tRdTpFAGZ3buVOX4qgYAgkl63CQmc0QdyYuOJU3G73VxyySU89dRT/PKXv4z8oBpNF0VKKCwUrF0rGDQIYqVNa6ewHGtrITn59EdKSvjtpz7CCWuIV199lby8PB5++GFuuOEGPvzwQ+bOnYvX62XNmjUMHjyY999/n+LiYi699FKKi4uZNGkS7733Xvu9ARpNB2btWli6VJCTA0YMKVKnsBzbkhtuuIGHH36Ye++9l+rqai699FIqKytxu92MGzeO6667juTkZHr16sWMGTMYNmwYAEeOHInyyDWa2GfbNpg3z6BfPyLqoNWexNhwYpPsOu9wcnIyR48ePakL+be//W2efvppduzYwYwZM+q3n0lXcI2mK3LgALz1lkGPHm3bXaelxJARG7uUlpYCqkomIyPjpMYT5557Llu2bKGkpKS+aw9AUlJSu49To+koVFaqlJ34eEhNjfZowqPFMQJSU1N5/PHH+fGPf8xDDz1U33gixJVXXklubm793yUlJbr5hEbTCIEAvPuuaiSRmRnt0TROp5hWJySorPpTCQSIOFrdFOnp6dx55531fw8dOpRnnnmGX//619i2TXx8PDNnNizD/dxzz3HXXXdFOnyNpssgJXzwgeCrr1Sz2limU4hjY2k4gYCMSByb4vXXX2fXrl2sWLGCiRMnApCTk8PNN9/Mnj17uPvuu5kwYQLp6el15wxwyy23kJOTc3Yn1mg6IWvWwOrVKjIdKyk7jdEpxLEt+da3vsW3vvWt07aHGk+88sorJ22Pi4sL25RCo+nq7NgBBQUG2dlE1Gc12mifo0ajaXPKy9XaL926xWZkOhxaHDUaTZsSCMDbbwtsW1DnfWoxJSVFzJ/ftmvHhIjKtNo0zVxgDrAW6AccsSzroWiMRaPRnBlFRUUUFhaSn5+P1+ttcl8p4cMPBXv2qNLAs6GkpIgnnmhYO+aOOxbW93RsC6Llc+wO/MOyrPcATNPcZJrmPMuy1kRpPI3y0UcfccsttzBjxgySkpIIBALcc889YfMY586dS35+PhkZGVEYqUbT9hQVFTF9+nT8fj8ej4eCgoImBXLDBli1SgVgzpbi4pPXjikuXtqm4hiVabVlWatDwnjCOMIk40Sfiy66iJycHH7605/y4IMPMnnyZH71q1+F3Xfu3LmUl5e38wg1mvajsLAQv9+Pbdv4/X4KCwsb3XffPlUa2KdP6wRgQmvHGIYLl8vDsGFts3ZMiKhHq03TvBZYaFnW5qb2O3jwYP3/HceJqCVYU+tLnwlSSoLBIIFAgEsvvZRf/epX/PrXvyYQCLBr1y5+8Ytf4Ha7Wb9+fX3HnmAwyIIFC3C73UycOJFrrrnmtOM6jnPSdTVGZxLcznItXfU6xowZU79uUlxcHGPGjAn7Ha6pgddfT8a2we+XtMZS7j16jOSWW96mpGQ5ubmT6NFjJMeOVeD3Q21tLQcPtq59FVVxNE3zIuAi4D+a2zcrK6v+/xUVFSctbNUUa9asidg/0hhCCNxud/05u3fvzpAhQ7juuuvYsGED//d//8djjz1GXl4et99+OwMHDmTlypU8+uijxMfHc8UVV3DdddeddlzDME66rqaIdL+OQGe5lq54HdOnT2fBggVN3lOOoypgHKf1E73z8qaRlzftpG0+H1RVCbKyGl9PviVETRxN07wSmAL8HOhjmmaOZVkrWvMcK1eu5KqrrorYPxIp5eXlxMXFcd9992HbNkePHj1tn+TkZB566CG6detGWVnZWZ9To4kVvF5vk/eRZcHGjbFfAdMcUfE5mqY5Dngd8AIfAe8Bw1r7PMuWLYvYPxIpCxcu5Gtf+xp33XUXDzzwAD/4wQ/qn3O5XEgp2bx5M/fffz/XX389//mf/6mbUGi6DHv2wKJFKtE7lnoztoSoWI51UemUtj7PlClT8Hg89ZZjfn7+GR+jsLCQXbt28cwzz5CcnMzx48f5zW9+gxCC22+/nYyMDD7//HNKSkq48MILeeyxxxgxYgTXXnstDz74IJMmTaK0tJQlS5Zw8cUXt8FVajSxQVUVvPOOIC2t4yR6N0XUAzJtyQUXXEBBQcFZ+Rzz8/MpLi4+bfsf/vCH+v/PmTMHgNzcXL73ve/Vb7/xxhsB9HIJmk6P48CCBYKqKujXL9qjaR06tThC8/4RjUZz9qxZA1980fH9jCfSwb0CGo0m2pSWwsKFncPPeCId8lKEEB1+bWgp5UnLLWg0HZGaGuVnTE0N72dsz1ro1qZDTqvj4+MpLy8nIyOjQwqMlJLy8nLiO4PXWtNlCTWurawUYf2MLamFLikporh4KcOGTW3T0sBI6JDi2Lt3b/bv309ZWVmTFqTjOBgxaOcLIYiPj6d3797RHopG02I++wzWrWvcz3imtdDt3ViiOTqkOLrdbvpFEBI7ePBgp6li0GhiiQMHYMECg759G6+bDtVC27Y/olro9m4s0RwdUhw1Gk308Plg7ly1cmBiYuP75eZ6ueOOhRFPk89UTNsaLY4ajSZipISlSwX790eWtpOb643Y+jtTMW1rtDhqNJqIKS6GFSsEAwa0zfHPREzbmtiLVmg0mpikvBz+9S9Bz57g7gJmlRZHjUbTLMEgvP++akOWmhrt0bQPWhw1Gk2zrFwpKCkR9OkT7ZG0H1ocNRpNk+zcCUuWqETvDlhz0WK0OGo0mkaprlZ+xvR08HiiPZr2RYujRqMJi+Oo8sCKCkG3btEeTfujxVGj0YRl/XpYv17Qt2+0RxIdtDhqNJrTOHhQtSFrrWVVOyJaHDUazUn4/ao8MC6u6fLAzo4WR41GU4+UUFgo2LdP0J49W2Kx72MXyHPXaDSRsmULfPqpoH//9jtnrLUqC6EtR41GA6jywLlzBZmZ7VseGK5VWSzQqcXRttVDo9E0zYnlgWlp7XvuUKsyw3DFRKuyEJ16Wn3kiEFFBcTFQUKCci57POrvuLiule2v0TSGlFBUpMoDBw1q//PHWquyEJ1aHG0bUlLUFMHnU4sBOU6DKCYmQlKSEk6Pp2t0GtFoTmXnTvj44+iWB8ZSq7IQXUIODOP0ldGkVFOJsjIlmKCsydRUJZjx8V03v0vTdaisVNPptLSuVx7YHJ1aHDdvdgOCgQMl3btDcnLDc0I0TK9D2DYcO6YEUwglkqmpysI8cT+NpjNg2/DRRwkcPy7Izo72aGKPTiuOK1cKVq2KRwjBtm0CjweSkyU5OZCVJcnIOFksQVmKJya9+v2wb58SysREyMhQgqktSk1nwLJg61Y3I0dGeySxSdTE0TTN3sAcIM+yrPGtffx773VTWJhR/3dKiiQzU9KtG/TsKendWzJ4sMPQodC7N6SlSdLTT/Y7ejwNU40ThTItDdLTla9So+mI7N4NH35o0KuXQwyuXhwTRNNynAy8B5zbFgfv2XMvGRkOtp1FTY2HqiqDqirBzp0n75eWdpyBA+MZNMghN1eSmyvJzoYePZSQhizJkFBKCVVVavodHw+Zmcqa1F8wTUehqgree0+QktL8LKikpCjmosjtRdTE0bKsN03T/FpbHLuoqIj586dTW+vHMFyAAHoAg4ARCDEWKccC51NRkcKGDbBhg1K3pKQyzjmnBtPsTU6OJCNDTcX79lV+S8NoEEy/H/buVdZmjx4qMq5FUhPL2DYsXCiorFTR6WPHGt83VitX2osO43M8ePBgxPvOnz8fvz+UcV8Xiqa07vEJUobyFQQwCsgHpgFTqanpzoYN3dmwAZKSfIwdKxk2zEdWloPHAwMGBOnfP0hmplOf9hAMwoEDApcLMjMdkpNlq4pkeXl56x0synSWa+mo17F2bRyrVsXTr5/DsWNw/HhNo/uuX7+ovnIlGPSzfv0ievSIPQel3w+1tbUcPFjdqsftMOKYdQZV8FdccQVPPfUUPp8fl8uFEIJgMIDjOBiGgculotjqg98IbASeRhUMeYGrgWupqRlCUREUFSXQo0c1U6cm4HI5HDwIHo9k5EjIyZGkpKjz2rbKpaypgZ49lSXZWnljZ3L9sU5nuZaOdh27dsGaNQZDhpyc2paeHr4kJi/vMj788Als24/L5SEv77JG940mPh9UVQmyspKb3/kM6DDieCZ4vV4KCgp4/fUF5OdfTlwcFBUV0q1bJkePHsHrzQfg7bdf5s03X8K2g/XT70DgU2A5cBdgIsS3kfI7HD7ci7fegvfeg4kTBfn5NuvXw7p1gv79JSNHSnr0UKk/waAK3iQmKpHUgRtNtKmsVHXTqamn5/w2RqxWrrQX0YxWTwW+C/QxTfMe4AnLso631vG9Xi+JieeQnNwDtxvGjPFiGMonGLLmzj/fy8yZN1JUVFgvmG+91SCY8BlSrkUJ5TeAnxAMfp1ly2DZMoPc3CN885vpHD4MCxYIevSQ5OVJ+vRRVqPPp6oPMjJo92J+jSaEbUNBgaCq6szzGWOxciUYhNJS2LHDoLQUdu0SBIOp3Hef0/yLz4BoBmSWAm3afiM1VZKYqCLMjqPe1Npa9ZyU6t+RI73k5Xnro3bnn+9l1iwlmCkpmTzyyJ34fLXAu3WPEQhxG1J+n5KSTH77WxgxwuHqqx18Pli8WInkeedJevVSEe7KSvXIylKWpa7p1rQXqm4aNm8W5OSc3bGiFbn2+WDrVsGWLeqxc6cgGDz5JnK5RF3gqPXO26ltmZQUGbZhZzDY8KiuVj7C43U2a1wcnHeel/PPVx/+6NGjTrEmtyLlz4D7gZ/jcv2CL79M4ssvDUaPdpg1y8bnEyxaJOjXT3LuuSrKbdtqql1VpabauuJG0x6UlMBHHxlkZ59dJkV7R66PHIF16ww2blSCeKoY9u4tycmR9O8vycqSuFyVuFxJrTqGTi2OjeF2N0xxQ8GUQED9QlVWKgGT8mShnDXrRlasKMQwMvn97+8kGCzDMB4CngJuA27n889T+eILwYUXSq65xuboUcG8eYJhwyRjx0pSU5XlumMH9OqlrUhN21JWppZVzcg4vW46ZAVmZ5vk5U1r9ljhei62VBwbs0CPHIFVqwzWrBHs3Nmg5EJIcnIchg+XDBsmGTxYnlTd5vPB3r2y1e+lZsVRCHG1lHJu65429gjVWaekqCl4ba3KAausVM+PHdtgTY4fP4oPPijk8893U1T0AnAP8HvgPqT8CcuXx7FyZYApU0r55jf7s22bYMcOgWlKBg6USNlgRWZlaV+kpvWprVXCGAgIevQ4+bmTrcA47rhjUbNCF+q5GIpct7Tn4qkW6O23f0BZ2QUsX26wZUuDIMbHS8aMkeTlOYwapQyL9iaS2/IeIcRk4EUp5edtPaBYwDBU1UtSkpoCV1erX+HaWiWgpunFNL1YVhHf/e7LdT7Jw8DtwP8CjxEMXsVHHw1mw4YafvADDwMHSj75RLB9O0yYIElLU9P5HTugb191Lo2mNXAc1YJs167wfsYTrcBgkIiswEgi15H4JBvOnUsg8DP++MdxBAJKhuLiVEDzggscRo6UUe8SFIk4fg/YCXxfCHErsEBK+a+2HVbs4HarOuq0NCWOZWXKmnS7Ydw4L6+8UsAbb7zM22+/RDAYRIgSpLwGKS8D/sCRI0N57DGYPNnhm9+0KS8XvP++YNw4yZAhEsdRda6ZmdRX4Gg0LUVK2LBBNV4ZMCC82+ZkKzAurBUYTuiailxH4pOUEjyeGcAFwKWAcmcNGuQwebLD+PEyplY7jEQcXYAN+IALgRwhxGXAMinlP9tycLFEqDNPdrbycRw9qqbdI0Z4eeQRL9dddyMff1zIsWOZvPbandj2IqQcA/w/4G4++SSejRsF3/2uzciRklWrBLt3g9erksjLylRQqHdvHazRtJzdu6GgQK033Zi75kQrMDvbPE3EWhJ8acon6Tgwb14JS5dmcOzYGGAMLleA0aMPM2NGj3ZdzOtMiEQcXwZSgfeBf5NSbgUQQjwEdBlxPJH4eCVi3bsrUauogFGjlE9SSvB6R/H7389h27YlSPkAQvyT1NT3OHbsHP70Jzder8O3v91gRU6c6JCT0xCsyc7W02zNmXPsGLz7riAxsfn1pkNW4LFjFac915LgSzifpG0rC/a99/yUlQ2v23M/U6ce59pr+5Gc3KPJY0abSMSxGPiRlLIytEEI4QG6tdmoOggejxLJbt1UpK2qSm2bPt1Lt2738P3vLycQ8GMY26iuHgv8BJhDUVESW7YIbrrJZvBgydKlBkOHSsaNk7jd6te/Z091XB3N1kSCz6cCMMePC/r2PbtjtST4cqI1OnToVI4encjf/+5i/36BkpntwG8Q4hW6d/8vkpPvOrtBtgORiOMNUsqT1vCTUvpR+SsalCXZt6+aFh88qERy/Hgvr75awJIlhXz22W5WrHgB+B0wj4SEdykrG8GTT7q47DKHGTMctm8XHDoEU6aovpKHDqkvfAcr39VEAduGpUsFX3119one0PKywdxcL37/RF57zcWuXepXvWdPyfjxW1i0aDyOUxNTqws2R7PieKowahonMREGDFABmwMHYPhwL+edp6LalvUyfr8P2EJt7bkYxn1I+V8sXOhi48ZqfvazRAIBwfz5gkmT1DS7uhr27NGdxzWNIyWsWwcrVgj692+9mcaZlg3u2wdvvOFi40YVUczIkHzjGw6TJjm43YMZO3Zeh6vR1rHRVibUKXzQIJXkXV0NeXnKirzggosRwgD8SHkfQlwE7GLv3jQefFAtjdmtGyxdarB2rSAhQTmz9+xxUdN4ZylNF2b7dli0yKC2togPPvgtJSVF7Xr+mhp4/XWDBx90s3GjQUKCKoCYMyfI1KlOfVAoN9fLFVfc1WGEEbpohUx74HarKpi0NNi/H4YO9XLnnffw3e8ux+/3AwIpP0E1Qn8Bn+8a/vxnuPRSm2uucdi0SVBWBhdeKPF4JLt3K/9menqUL0wTMxw8qDp6HzlSxLPPtm9TWimVtfrmmy4qKwVCSPLzbWbMcEiLva5mLUKLYxuTmAg5OSqqPXSol7/9rYA1awqBTP7whzsJBisQ4lsMGvQSX331TT74wMVXXwluuUWVHy5YAGPHGmRlKZH1+1XXcR2o6dpUVakAjG0L9u5tvdK+SCgthVdecbF1q5p4DhnicP31NgMGtNkpo4IWx3bAMBqWUYiP9zJ6tJekJMjLG8Vzz73C8uUvsn37DRjGn0lK+hfbtqXw4IMBfvYzF1lZsGRJAvHx0K+fEtlAQFmROmG8a+LzqRZ5Bw4oP2NrlfY1RyAA8+YZLFhgYNuC1FTJddfZeL2tX9ccC2hxbEcSElTApqxMpf5MmODls88K+eSTIFLaOE4htbXDgNeoqsrn0UcdvvMdhxEjHD76yGD8eMnw4bI+UNO3r67L7mrYNixbJvjyS1FvqbVHU9pt2wT/93+h1BzIz7eZOdM5bXnjzoS+tdqZkBWZnKwW5zrvvHw8Hs8Jfsj9qPVsnsBxbuellwwmTUri29+GVasEFRVgmhK/X7W979fv9I4rms6J46i1plesEKe1IGurprQ+H7zzjsGSJQZSCnr3lsyebTNkiGz1c8UaWhyjRGIiDBwISUlennuugM8+KyQ5WTXXDQb9GMZ/cs456Wzb9l2WL0/kyBGHm2+22bZNUFUFkyerL+fOnUogY6kmVdP6SAlbtsDixUa7lZhu2yZ44QUXhw4JDEMyfbrNN77hdJnyVi2OUcTlUr7Dyy5T3cg9HtVc97nnXmHJkhfZsuWHGMbfiI+fx+bNycyZU8svfxnHoUOCDz6Ar31NEh+vKmr69CEqbZ007cOuXSoAk57e9msSBQLw3nsGixYpazE7W3LTTcFWSTDvSGiXfpQRQqXn5OSoaZNatqE/Ugbroo/L8PvHAJ9RVpbEnDlQViaorRUsWCCorFRW4969Ta9BrOm4HDyoFsfyeESbp8ns2QP/8z9uFi5UlQdXXGHz3/8dXWEsKSli/vz2z+HUlmOMEB+vBPLAARg7VvkhAwE/UgocZycwGXiV2tqr+d3vlN9n7Fi1QPtFFzn1qT6BgGp/1hmjh12R8nKVy1hbK+jdu+3O4zjw4YcG77xjEAwKsrIkP/yhqv2PJu29PMOJaMsxhnC51PT40ku9PPtsAT//+f3cd98fiIvzIEQthvEtsrM/xLYFf/ubm8WLDVJS1Jd6926VKnTkiLI0ZOf3l3d6qqvh/fcFR460rTCWl8Mf/uDijTdcBIOCqVNt7rsvGHVhhPAdgtoLbTnGGEKobjzf+Iaqy3a7ITu7H//85/ssWfIie/dejmHchpRP8v77Lo4cEVx/vc3SpQYTJqg1NsrLVcqHzoXsuIRyGXfsEG3a73DTpjj++U83VVWClBTJ979vk5cXfVEM0V45nOHQ4hijJCWpaPbevTBkiEle3noWL1Z+SCGeYuTIwWzefAsrVngoL4cf/9hm5UpBbS2MHatyIffuVZaoblzRsfD7YckSlcvYr1/b/MAFAvDWWwaLFysn5siRDjfdZJOR0frnOhvaI4ezMbQ4xjBxcSpNp6JC1vsh/X4/QrgoLv5PbPtlYC5fftmLJ58U/OxnQTZsUAI5YYLE52tIFu8q6RcdnWAQli8Hy1K5jG3xw3bwIPzlL2527VIpOjNnOlx6qROzs4y2yuFsDi2OMY7LBVlZDtnZyg+5fn0hO3fu5q23XgBWARfi8XzMnj39efRRN7fdFmTrVoHfr5pWBIMq1Ucni8cORUVFFBYWkp+fj9fbcNMHg7B6NSxfbrTqD9qJ68GUlU3kxRdd1NYKevSQ3HBDBWPG6Lbz4dDi2AEQQi3JMH26l3PP9bJpUxHvv/8yPp8PKbfj949DiLmUlXl57DE3P/2pzZ498NFHkJ+v/Eehapq2zpHTNE1RURHTp0/H7/fj8XgoKCjA6/Vi2/DZZ7BkiUGvXq33QxaK9gYCAsPoheNMAWDcOIfZs20CgWDrnKgTEqOGtCYcqamqNnv0aC8vvFDApEmh/pCHgEvp06eYmhrB73/vYu9eFeVcvFhg29Qvv6D7QkaXwsJC/H4/tm3j9/spLCzEtuHzz+GDDwx69GjdH7Di4qUEAv2BZTjOjzGMIDfcYHPLLbZep6gZoiaOpmleYprm06ZpPmCa5v3RGkdHIyFB5UOOG+flJz+5h/j4eCbTj7oAAB8DSURBVAzDhctlk5//EQMHriUQEPz5zy42bTKorhYsWiTw+aivpqmsbP48mrYhP1/5jl0uFx6PhylT8vniC1iwwCAjoy0WVpuJcr+cD2xn9uwNXHSRo/NgIyAq02rTNJOAZ4BRlmX5TNN8yzTNaZZlLY7GeDoaoUBNXJyqy16/vpBgMJM//emXdXXZD+A49/Dyyy4mT97DjBl9WbgQpk2TpKWpKLZunBsdvF4vBQUFFBYWMmVKPmlpXubNE6SmqjzV1iIYhLffNvjgg5EA9OmzmW99q4xRo8Y3+ppwa1V3ZaLlc5wI7LQsy1f393LgSkCLY4SEEsbz81Vd9osvPorj+AEbKe9HiL1I+RSffJJDdfU+brihBwsXCi65xKF7d1VNY9t6hcNo4PV6mTDBy+bNKsk7OVm0al380aPw7LMutm0zcLkks2Y5XHJJLkLkNvqaaFaixCrREscs4MTJXUXdtkY5ePDgGZ+kvLz8jF8TizR3HVIKcnLOJS4ujkAAVOuzvwB7gX+wbl0fgkEfM2dW8/bbgkmTfGRl2RQXC9LTJZmZ7TfN6iqfSVM4DpSUuPjwwwQ8HpBStlpd/Natbl55JZWqKoPk5FrGjHmRrKwBVFSEtxiPH1dO6PXrF9VXogSDftavX0SPHiNbZ1BtzLZtq/n886UUFFzAuHHjWu240RLHg8CJv5VpddsaJauFa5S29HWxRlPX0asX9OlzGenpC/jss0JSUkKtz95HiCswjH+xcWMyPl8cP/yhzbp1KUyZ4jBggGq37zjqGO2V59YVPpPGsG3YvFn1ZMzMFGFdGy2Z3jqO8lu++65R92N5lNLSc1m5spQ1a5q2BA8f3kR19QFcLjeOAy6Xh7y8y0hPj/3FYEpKivjrX2cSDPpZs6Yh+t8aREscVwA5pmnG102tJwFPR2ksnYK0tIaSw7g4GD58FM88o5ZgCAYvBArYsqUvTz0luPXWIEuXGni9kqFDJZWV6qbV1TRti23Dl1+qssCEhMaF8UyntzU18MILLtavV79uV1xh4/E8za5dpc2uK7Njx2r+8peZdb5qF1Om/JCJE2/sMFPqE2uvQ9H/1hLHqESrLcuqAW4F/mia5hxggw7GnD1JSSrVx7bhvPO8mKZqfQYbEGIyiYn72bNH8Nvf2kgJRUWCjRsFycnUV9OoabmmtQkGVbpOU8IIZ95oYdcumDPHzfr1BklJkttuC3LttQ7Dh+fjdnvqMhkar0kuKVlefz4pbbp3739GwhitdmIhQrXXQqjof35+fqsdO2pJ4JZlfQB8EK3zd1bi45VA7t0L48adWHJYit8/DnibiooLePTRALfdJli7Vgnj+edLAoGGZPH4+GhfSechEICNG1X3pORkmuzJGEmjhdC0u7b23/jwwyEEAoIBAyQ/+UmQiooi5s9XU/JIapJzcye1uLFDWwdxInEv5OZ6ue22haxdu4gHHrik1axG0BUynZJQqo/L1ZDqs2nTbgoKXgAuBt7g+PErePJJm1tuUVO92lrqV5HbtQuys9si567r4fOpypePPzZIS2s+Xae5RgslJUU8/vgMgsHHgKEATJni8O1v2+zadbpYXXHFXU2eb+DA8S1u7BDOym0tcTwT4R082EtCwii83tZd7UuLYyfF5VINJyZP9jJ2rJctW4pYvPhlfL5a4GrgWYLBH/D005Ibb7QxDHUjT57csPSCzoU8O44fVwtiLV+uErzDrdQXzjpqqtHC2rXrCQY/RCV1H+f88z9k9uzLgcjE6kzP1xRt2U6sLYU3UrQ4dmIMQwmc2w1Dh3p5+eUCHn98DqtWLUHKHwL7kPK/eeklN5Mm7eSKK/qyeDFMnSpJTtadxc+GqiooKoLVqxsvCTzTaelnnwmWLr0FiAO24XbfyGWXPVb/fHNiFe58Z5Ou05iVG2m0van9otnHMYQWx06OENCzZygK7eWOO+7hu99djt/vxzAeQsp9OM4fWL48h4qK/dxwQyYLF8LFF0vS09Ua24FA+6b6dHTKymD5chXsaqqJRKTWkW3Du+8aLFigUgmGDDnM0KHvMGbMYyft39yUPNz5zjaX8VSrM1LBb26/aPZxDKHFsYvQvbuyIMHLiy8W8OmnhaxZs5vly58BdgCvs3Fjb156yeE737FZsEBw8cVqbZqqKiWQffrovpBNIaXqlfjxx4KSEtHsEqqRWEdHj8Jzz7nYutWo77142WXpCPEfYY/Z1BS5PayxSAU/kv2i1ccxhBbHLkRamhJIIVTJ4WefFbFq1csEAgsQ4jJcrgI2bUrj6acFP/5xkEWLDKZMccjJUQGb0BrZuu3Z6di2SoVaulRQWiro27f5nNHmrKNNmwTPP++islJVMt18s83QoS1fwiDc+Y4dq2jx8cIRqQDHwrS5ObQ4djFCuZB79kBenvJDhtbJDgTOB96ntHQ4v/udm5/8RCWLjxsnGTlSNc7dtUv5Mdt6idCOhN8P27ZBYaGgqkp18I7URxvOOnIc+Ne/DObNU9UuI0Y4/OhHdqu8521tjUU6HY6FaXNzaHHsgiQkKIEsLYVRo7zk5RWyeHEQKEGIyWRmWhw+PJDf/c7NTTfZWJZqczZ+vCQpCfbtU5Zkjx7aD1lVBZs2qXJAKZWP8WwoL4fnn3dRXGwghOTqq22uvDJ2lzAIR6QCHO1pc3NoceyieDzQv78Sury8hmRxl6uG66/fw7p1A1i+3ODZZ13MmOFgGA6VlSrVJyUFjh1TAtlV/ZBSqmVw169XEemUlOat6eaiuBs3Cl54wUVVlSAtTfKjH9mMGBE7KwF2NbQ4dmHUsq8qF/LZZwvYsKGQESPyWbvWISXl10yd+h0KCwfy7ruqs/hVV9kUFMBFF0m6dVPiuGOHyqcMl8PXWQkGobTUYMcOQXGxoGfP5iuKmorOBgLwzjsGH3ygnJQjRzr88IetM43WtBwtjl2cUC7kRRepIE1xcRF//evl+Hx+DONhRoz4HVu33sKqVS4OHIDZs20KCgSTJ6uuPsGg8l92767yITvS9K8l1NRASQksXRqPzycibtbRWHR23z547jk3u3cLXC7JjBkOX/96x5pGd1b0R6Cpz4XMyoJPPlFrnIRu4k2bfoaUXpKTK9i50+CJJySHDws+/tjgs88EhqFK4o4eVVU1Pl/z5+uIOA4cOgTr1qk1pauqjDNaOjUUnQ01ghg6dCoff2wwZ44Sxp49JXfdZTN9uhbGWEFbjpp6unWD6dPzefZZDz5fLVJKQBIMrsW2hwCvUlMzjT/+0WHWLAk4HDkCEycqP6TPp9J9evVS/rfOUlVz/LgKXm3aJNi0SZCRAd26nZkv8MTobL9+lzB//gQ2blQq6PU6TJxYxJdfLsFxwleaAC2uRNG0DC2OmpO46CIvc+cW8PzzL/Puuy9h20FA4DiHga8Dv8Vx7uCNN2D7duWHnDdPMGWKU5/0fOAAVFTQqkuMRgPHUUGXnTth/XrB4cMqGu1206LO3bm5Xo4encjf/qaCLklJqq69e/cVp/kjgfpthuFCCIFtB8M+r5c1aBu0OGpOIz/fy8SJXmbNupFPPy0kKyuThx66k0DAjxD/j6FDkykp+RFr1qgp4fe+pxLG8/Iko0c3WJE7dqh0n4yMjueLrK5Wbd9KSlQZYHy8Cjy1lMpKeO01F6tXqzdixAiH73/fpnt3mD8/fA/H0DbHcZASQIZ9PlqNGTo7Whw1YYmLg6uuUgtBVVbC0KGjePPNl3nrrZfYvPlnuFxPk5HxEQcPZvLkky5mzpQYhkNpKVx4oSQjQx3j8GGVu9erV8eIaPt8ylosLYUvvhDs36+i0S21gKUEyxK89pqqdPF4JNdd5zB1asO6PY1Vi4S2CaEsR8cJhn0+VitMOjpaHDWNEopkx8eDlF769SvEcYJ1izB9wbFjA4FHse1beeMN+PJLwcyZNu+/Lzj/fMmwYcqKDARURDs5WU1VY5FAQDWMOHwYdu5UvkWPR6U6Nefba+z5sjJ4/vljbN3aA4Bhwxxmz7Y5demZxqpFTtwGp/scY73CpKOjxVHTJEKoNJ24ODj3XJUsHgj4UX7IGuDfEWIJbveLfP55Ijt2CK6/3gZg+3bVQLdHD/V6nw/27XMTF6eOGQvdxgMBFWkvL1cJ8Rs3Co4fF2RmqjGfmp94/fVPUFV1hOxsk7y8aWHzFwcO9PLxxwZvvy3x+3sAFRjGf3PNNf9GVlbj5XTNNV6ItcYMnR0tjpqISE2Fq6/24nIVsHKl8kM+/PCddUswzOW88x7g0KG7+eqrbjz/vJtx4xyuuspm/nzB0KGSsWNV6WFKisPx48ofmZKiIuSJie0f2a6tVaJYWan+/fJLNYUuLy9i584Ga+zE/MRg0Merr96O40jc7jjuuGPRafmLn35awssvT2bPntAFvQP8DDjAli29OOccLWYdBS2OmoiJj4cZM7x4vcoPOWTIKN5992XeeOMlVq9+ApfrKS67bANLl+ayZo3B5s2CWbNs3G4V2R4zRrJ37xpKStbj9eYzapSXPXtU9DczUzXFaMtSxGBQBVrKy5UVW14OW7cK9uxRi4zV1hbxv/97shV4oj9QWcsOUjoEgw3TXLfbQzDYEykfobDwBgAyMyVTp27iX//6jvYLdlC0OGrOCJdL1VMnJgJ46d27ENsO1kVV/TjOP7j33rt4/vkqduzoxosvujnnHIfrr7eZN28lL7xwPY7jx+Px8PLLBZx/vpdgUPVBlFIdNy1NNcfweM7OopRSdcyprVWpRbW1SiCPHYPNmwUHDggSE9X1CAHLl58eNb7iirvqfXspKZn84x931IldHMOGTSU728uECV9SVNQX247DMIJ4vXu54YY+xMcPZejQ2PcLBoMND8dR71vINyxEw8Ptbnh0tOyDlqDFUXPGCKGmwwkJDSscBgJ+4uI8TJuWT0lJEXv2fB34JvAY27b1Ys4c6NPHIRhMAw7g8/l5881l5OR46d69IZLt9yuhDJGUpB4ejxJml0vdmCeKZuhmtm318PmUEFZXN+zj98P+/YJNm+D4cWUpnpqa01jU+ETfXnb2aIqLl9Kr1wS++upC/vQng6qqnLr35U0c525Wr95Lfv7C+tfFgig6jnpfjh9X/9q2et8qKkLr26i1g5KTG95jw2h4X/1+9aishJoaUZdapD4Hl+vkz6mzoMVR02ISE9U0Oz6+gMLCQiZPzscwJO+++zC27QNeAt5HiDlIeTP79k0FihHiCVyuZ8jImMqCBQYpKZJhw6BXL5UCFLrBpFQBk8OHqb8ZQd2QJ96cJyKluqldrtBrBVu3wpEjqtQxI0M9whFJj8H+/b3s2HEhr70Gx4656l7nkJ39IsuW/QiwsW1X1PMObVv9ONTUNFiEGRmS3r1VqWi3bpLu3cHnq6Zfv6STBDFkKYaQsuGhxFJy/LgSyiNHoKxMsG+f+lGrqVHvc1ycmgEkJXXcSiktjpqzwu2GK6/0Mnmylw8+KOLHP56Oz+dDSgchDISoRMrbgD8BTwDTkfIh4uLu5fBhg5EjHYRQ62fv3LmS3buX4vXmM2WKl7Q0WW+RNHWDOY6yFGtqoKJC3aj790NtrUAIZQ316RPZ9TRm6VVXw9KlBkuWGBw7pgaTnS259lqbsWMl27fnsmJFg9WZkpLJ/Pm/bdfpdMh94Per9ysrSzJkCPTtK+nTR70PcXEn14MfPCibXYI33HufmKgyDnJyACRSKhGuqJDs2wd79gh27IBdu9RnkJDAST98HQEtjpqzRgj1xd+6NdS0wsEwDCZNupjJk6/l8cfvJBjcgts9iyuv/JiNG01KSuJ4802YN8/ga19zGDzY4u9/V8GQxYs93HTTQgYM8NZbgomJkoSEUM6lOm/DNFGcZFkmJChBPXKkIf8wJSXy/MQTUUsfGHz6qYHfr1RiwADJRRdVceGFCfW+txOtzpBvsj1K+6qrlQUXCKip8ZAhanwDBigxPFu/baQIoYQ3M1M9Ro+W2DZUVkp274biYsH27eqzio+PnVSuptDiqGk1Lr44n8ceU01z3W4Pt956D16vl/POG8W8ecuIj59KVtY5TJ9u88UXDgUFBlu3GhQUuBBiPFK+DPyZYHAphw4tZfx4JSjKlygIBk/u+uNyKQsmNfX0AEFzq9s19Xx1NaxZI/j0U4OSkoYDjxzpcNllDiNHSioq/BjGyYvphKzO+fN/26alfbW1KqgUCEBqqmT0aMjNlWRnqx+GSDsFtTUuV4MbY8wYSW0t7N0rKS5Wvt+aGkFCghLKWLQotThqWg2v10tBgfI/er355OR4qaqC88/3Yprq/0uWVLJvn7qZf/Urm5ISh0WLDNatE8BMYCZS7mPPHptt2wSDB8v6QEw4QtZfSkomVVVH6oMoc+c+RDCopvfhBOrU/MQNG1axf/9E1q0z+OILQTCozK3ERInXq8r9srMjex/aYvGoQEClHvn9EB8vGTEChgyJPUFsioQEGDwYBg+WTJsGu3dLvvhCCWUgIEhLi606fC2OmlbF61V5kKB8gceOqYCKy6WSvi+80Iff77BqlWDfPkH//pJbb7UpL4f33tvLunVpVFf3YfVqWL1aWUajRkmGDnVITFzHgQMLGT5cTYND1l8g4AOUj9Mw3AghCAYD9dvCCdSAAdMwjNXY9iSkvJiCgnFIqQRRCMnIkQ4XXOAwbpw84+lfay0epaLJatrscsGgQZJzzpEMHqymzO4OfPd6PJCbq34kL70USkoka9cKdu4U9XmvKl0serT722uapgH8GHgYuNiyrM/bewya9sEwVMpPcrISyKoq5bTv2xf69VvB/PnL2Lt3KlLCgQNLmTx5KrNne9m+PciHHx5g06YkKiu7UVQkKCoygPFALkJs4PzzD+L3HycQuBzYC5Qh5TFs2wFsIA5IY+DASxg37qeUlg5n3TpVBbN3r+DQoYnAe4ASIZdLMny4w7nnSs491yEjQ1mlixe3TODOJoXH71dVO4GACqqMHSsZOlS9l7Hup2sJiYkwejSMGiU5fFhZk2vXwoEDgvR0SE+PjjUZjd+ePGAlUBOFc2uigMejBLG6WjWMXb5cRbWVb9KFlIJgMMi8eR5+8YuFGAasX/91AgE/bnceU6f+ky+/9LN3byaQhZRfY80agGzgsibP/dVX6nH6mCQDBkgGD5aMGCEZMuRkC7E5n2Vr4zjKSqyuhrg4FVgZOlTSv78KLsXKVLMtCXWk/9rXJBdeCNu2SVatUtZkfLx6rj0Xc2t3cbQsax2AaZrtfWpNlElOhv79bebPV1Ftx7EJBFQphpSqV+HatUtxuVSvQpUz+Bn79v2EadNm8tprd2DbPTAMk/z8x4Ecdu0qp6yslmAwGZ8vDsPwIKVASh+JiQapqQmkp0sM4yCBQDE9e/qIi/sK0xzJkCGNi11ja760NicGV/r0kZx3nrIS09M7p5UYKR4PjBwJI0ZI9u9XU+4NGwSOo3qENpd+1Bq0iTiaprkQCLeC732WZc1tyTEPnlg2ESHl5eUtOVXM0VmuA6CiopwJE8bg8cQRCIBhqEiC49i43XFcffVYdu92s3ixpz6g8uWXS9i69RNmzPg1NTVl5OZOYuDAdCD0viSgptI2UHvKGWvYsWM1f/nLtQQCfjZvVn7IZcs83HLLOwwcOD7sOLOzTdzuOIJBcLniyM42OXasov754uJl7Nmzpm4s4Y/RGLat8jEDAdXfcfDgADk5Nn36OCQlqZyklnQabwkd4bvlcsH48TB8OBQXx7F2bTxVVcofnZ4u60pEazl4sLr5g50BbSKOlmV9vbWPmXVqE7w2fl2s0VmuA2D69OksWLCAwsJCJk3Kr4tiF2Ka+Uyc6MUwYMyYAp54Yg4bNiypizgHcJxqZs6874zPV1pq1QdogPrjlZZa5OVNC/uavLxp9V13TvU5lpQU8fe/f5tgMBDxlFtKFVgJlTSec45k0CDJoEEqQhtNK7GjfLeyslQQ55JLYNs2WLFCUFoqiIuDhARBVlbrdlPuwPEuTUfmxKg2wLRpXioqVDmalOr5e++9hxtvXI7fr1JiMjOnUlEBBw8WsWVL5IGShs45DZU7kaTYNBZUUVPuQLNTbimVGFZXK2uxRw/J8OGSgQPVjZ6U1DFScGKN0JR7+HBJaamkqEiwf3/rd1GORrS6G/BTIB242TTNVy3LKmrvcWhiC7dbJQOnp6uo9pEjMGSIlxdeKGDt2kImTMinT58JzJu3gieeuBzbjjxQcmr1SigfsqU+RCW2cdg2p4ms46jx19Q01DOPHq1KDbOzVTpTLCY8d0QMA/r3h/79Jfv2HQdSW/X40QjIHAXm1D00mpNwuZRApqWpYEVampfRo5WIxcdDMLgUx2kIlFjWUlJSvCQnNx3Vbc3uOLm5Xm655R1KSy2GDZtK//5ejhxRQRWA7t1VtLl3b9XoITVVJUB31AYMHYG2sMD1tFoTkwih8t8SE1UKR02NClLk5Z3cIm327Cn06iXZsQMOHRL1tdhJSeq1rX3ThFqipaVNoFcv5a+sqVFWYa9eksxMFU0NCWJXSMHprGhx1MQ8breyJNPSoHdvL927F/DRR4Xk5eWTl+fFMJQfT0pJeblqU1ZaqtqU2XZDi7O4OPU4sWFrSLyEUNPgULPXQEAlrAcCDb0PQe2fni7JyQmSk+PQrZuydFNTVapSfLy2EDsLWhw1HYq4OBW8mTbNW9/AtaZGBT38fuqm1yoKLIRK8zh+HKqrBRUVyh9YVaVeFwiomt4ThS8uTuJ2N1ieiYmQkiLxeBoausbHw/HjfgYMUP/vyGV8msbRH6umw6JamalHZqay+pTgKX+l36/2U3Xdkl69Tm6aqzh5w4lWn8ulxC8+Xk2R3W4VTDEM1QexI6zDrWk5Whw1nQbDaBCzlJSG7VKqqfGJ0+aQSEqpBHH16iKWLy9kyhSVa+ly6elxV0eLo6bTE1ocqjGKioq45hpV6+3xeCgoKDgpB1PTNdGxNE2Xp7BQ1Xrbto3f76ewsDDaQ9LEAFocNV2e/HyVHuRyufB4POTn50d7SJoYQE+rNV2eEzuY5+fn6ym1BtDiqNEAp9d6azR6Wq3RaDRh0OKo0Wg0YdDiqNFoNGHQ4qjRaDRh0OKo0Wg0Yegw0er7778/2kPQaDRdCCFPr8TXaDSaLo+eVms0Gk0YtDhqNBpNGLQ4ajQaTRi0OGo0Gk0YOky0uilM07wEmAkcBKRlWQ+e8nwC8DhQCgwBfmNZ1pZ2H2gzRHAddwG9gf3AOOA+y7I2t/tAm6G56zhhv+8ALwOplmVVteMQIyKCz0MAt9X9ORDIsCzrB+06yAiI4DoGoe6P1cC5wKuWZc1t94E2g2mavVGrluZZljU+zPMG8GugCsgB/no2yz53eMvRNM0k4BngF5ZlPQCMNU1z2im7/Qewy7KsR4DfAX9t31E2T4TXkQL80rKs3wJvAY+17yibJ8LrwDTNEcDIdh5exER4HTcC5ZZl/dGyrF8Cv2/nYTZLhNfxK+ATy7J+A/wWeKJ9Rxkxk4H3gMZ6tP8bkGZZ1hzgLuBF0zRbvP5khxdHYCKw07IsX93fy4ErT9nnSmAFgGVZG4E80zTT2m+IEdHsdViWda9lWaHcKwP1CxlrNHsddTfsr4CwFmWMEMn36jtAd9M0bzdNM2SxxBqRXMcBoGfd/3sCa9ppbGeEZVlvApVN7HLifV4G1AKjWnq+ziCOWZz8hlXUbTvTfaJNxGM0TdMDfA+4px3GdaZEch3/AzxsWZa/3UZ15kRyHTkoS+WPwN+BBWdjqbQRkVzHk8AFpmk+CdwH/K2dxtbatOp93hnE8SCQesLfaXXbznSfaBPRGOuE8c/Af1uWVdJOYzsTmrwO0zT7A92AfzNN8+66zb80TdNsvyFGRCSfRwWwEqDOh50G9G+X0UVOJNfxd+D5OtfAtcDrpml2b5/htSqtep93BnFcAeSYphlf9/ckYJ5pmt1PmDrPQ00vME1zDLDesqyK9h9qkzR7HaZpJgJ/AZ60LGuNaZqzojTWpmjyOizL2m1Z1vcty/pNnY8L1PVY0Rluo0TyvVoMDAao2+ZCBctiiUiuoz+wr+7/RwGHDqINpmkmm6YZcgmceJ93BxKAL1p67E5RPmia5qXAN4FDQMCyrAdN03wUKLMs6zd1ovI46gtwDvDrGI1WN3cdbwOjgb11L0kOF7WLNs1dR90+PYFbgIfrHn+xLKs0WmMORwSfRzrwKLATyAXesixrfvRGHJ4IrmMyKmi5FhgErLEs65nojTg8pmlOBWYDl6NmT08APwDGWJb1k7po9SNADTAAeO5sotWdQhw1Go2mtekQprNGo9G0N1ocNRqNJgxaHDUajSYMWhw1Go0mDFocNRqNJgxaHDUajSYMWhw1Go0mDFocNZ0CIcR/CSGOCyGmCiF+KYSYL4QYGu1xaTouOglc02kQQtyLqozwA3dKKY9HeUiaDoy2HDWdif8BpgCbtDBqzhZtOWo6DUKIWaiuLL8CviGl3B7lIWk6MNpy1HQKhBA/AO4GClDNWt8RQkyN7qg0HRltOWo0Gk0YtOWo0Wg0YdDiqNFoNGHQ4qjRaDRh0OKo0Wg0YdDiqNFoNGHQ4qjRaDRh0OKo0Wg0YdDiqNFoNGH4/+obGuw3bkwDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "train_and_test_approximate_gp(gpytorch.mlls.PredictiveLogLikelihood)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
